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I.        INTRODUCTION 


A.   BACKGROUND 

A  Magnetchydrodynamic  (MHD)  generator  is  a  device  for 
converting  the  energy  of  a  flowing  gas,  or  of  a  liquid 
metal,  directly  into  electrical  energy.  As  with 
conventional  electrical  generators,  the  conducting  medium 
crosses  a  magnetic  field;  however,  MHD  deals  with  conducting 
fluids  instead  of  solid  moving  parts.  Anticipated 
efficiencies  are  far  above  those  of  conventional  generators 
because  of  the  higher  temperatures  in  the  conversion 
channel.  The  thermal  energy  of  flowing  gases  is  changed 
directly  to  electrical  energy,  eliminating  the  mechanical 
energy  step  cf  the  conventional  generator  which  means  that 
there  are  fewer  parts  to  wear  out. 

Although  MHD  power  generation  has  been  studied  for 
nearly  30  years,  financing  of  MUD  research  waned  in  the 
1960's.  Interest  is  now  being  rekindled  because  of  a  new 
"energy  consciousness"  of  government,  industry,  and  the 
military.  MHD  offers  an  attractive  alternative  energy 
transformation  means  which  is  relevant  to  the  pursuit  of 
clean  ways  of  using  the  vast  coal  reserves  in  the  United 
States  and  flaking  more  efficient  use  of  dwindling  petroleum 
resources.  MHD  promises  to  fulfill  the  reguirements  for  a 
light-weight,  high-power  source  for  military  use.  For 
example,  the  Navy,  recognizing  its  energy  reguirements  and 
dependency  on  world  fossil-fuel  supplies[ 1  ],  is  sponsoring 
research  on  practical  MHD  devices  with  large  power  outputs 
and  high  eff iciencies[ 2 ].  The  Air  Force  is  interested  in 
power  sources  for  airborn  weapons  systems. 

The  basic  MHD  device  consists  of  a  channel  through  which 
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ionized  gas  or  liquid  metal  flows.  A  magnetic  field  is 
imposed  perpendicular  to  the  flow  direction;  as  a  result,  an 
electric  field  is  induced  and  power  is  tapped  by  means  of 
electrodes.  Principal  concern  is  with  ionized  gas 
(hereafter  refered  to  as  plasma)  devices  since  with  these 
the  combustion  energy  can  be  more  directly  changed  to 
electrical  power. 

A  major  problem  involved  in  the  performance  of  HHD 
generators  is  the  high  voltage  loss  in  the  vicinity  of  the 
electrodes.  Accurate  performance  predictions  in  the  design 
of  MHD  devices  require  a  realistic  determination  of  these 
voltage  losses.  This  research  investigates  the  nature  of 
electrode  voltage  drops  and  computer  models  are  developed  to 
describe  the  physical  phenomena  involved.  With  the  results 
from  this  investigation,  appropriate  steps  may  be  taken  to 
minimize  the  losses. 

The  principal  voltage  loss  mechanisms  in  the  MHD 
generator  can  be  divided  into  two  main  classes,  ohmic  and 
sheath  losses.  Ohmic  drops  are  those  that  occur  because  of 
the  finite  conductivity  of  a  real  plasma.  Thermal  boundary 
layers,  degree  and  kinetics  of  ionization,  and  Joule  heating 
are  factors  affecting  the  ohmic  resistivity  of  the  plasma. 
Sheath  drops  occur  as  the  result  of  the  Debye  shielding 
which  forms  a  non-neutral  layer  adjacent  to  the  electrode 
and  results  in  a  space  charge  field.  Material  problems 
restrict  the  temperatures  at  which  the  electrodes  can 
operate.  In  many  cases  cooling  of  the  electrodes  is 
required  since  the  plasma,  in  order  to  maintain  a  high 
ionization,  must  be  hotter  than  the  working  temperatures  of 
most  materials.  This  temperature  difference  between  che 
electrodes  and  the  plasma  further  aggravates  the  voltage 
losses  because  of  the  presence  of  the  thermal  boundary 
layers.  As  will  be  shown,  voltages  losses  can  be  as  much  as 
50%   or   more   of   the  total  power  output,  with  sheath  drops 
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accounting  for  a  non-negligible  fraction   of   the   drop   and 
boundary  layer  losses  the  rest. 

B.   REVIEW  Of  PREVIOUS  WORK 

Analyses  and  experiments  on  boundary  layer  phenomena  in 
MHD  generators  are  numerous.  It  is  not  the  purpose  of  this 
work  to  further  add  to  the  literature  on  the  boundary  layer 
losses,  but,  as  will  be  explained  in  Part  C  of  this  section, 
this  phase  of  the  problem  was  treated  as  a  necessary 
appendage  to  the  problem  of  the  sheath  drop.  Conseguently  a 
short  review  cf  boundary  layer  work  is  presented  here:  High 
and  Felderman[ 3  ]  as  well  as  Argyropoulos,  Deraetriades,  and 
Lackner[4]  studied  the  problem  of  the  turbulent  boundary 
layer  and  several  descriptions  of  the  nature  of  the 
temperature  field.  Doss,  Dwyer  and  Hoffman[5]  investigated 
the  boundary  layer  and  used  a  rudimentary  collisionless 
sheath  as  an  "inside"  boundary  condition.  Kessler  and 
Eustis[ 6  ]  reported  on  experiments  with  electrode  temperature 
effects  on  turbulent  boundary  layers,  and  Rubin  and 
Eustis[ 7  ]  extended  Kessler* s  work  to  include  the  effects  of 
electrode  size.  Wu  et  al.[3]  and  Oliver  and  Mitchner[9] 
investigated  non- uniformities  of  current  distributions 
within  the  boundary  layer. 

Although  the  existence  of  the  sheath  is  well  understood, 
its  effects  have  been  investigated  to  a  much  lesser  extent 
than  those  of  the  boundary  layer.  This  is  principally 
because  this  loss  is  described  by  a  relatively  complicated 
set  of  coupled,  non-linear,  partial  differential  equations 
that  present  considerable  difficulty  for  numerical 
solutions.  There  is  no  known  exact  solution  for  these 
eguations. 

Most  work  on  sheath  phenomena  is  embodied  in  "probe" 
theory   ii  estigations,   that   is  in  the  mutual  effects  on  a 
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quiescent  plasma  which  is  disturbed  by  the  presence  of  an 
electric  prcbe.  Such  work  is  relevant  to  MHD  electrodes 
since  the  anode  is  essentially  a  heavily  biased  probe  in 
contact  with  a  plasma,  which  is  quiescent  within  the  sheath 
region  (see  Part  C  of  this  section) .  The  cathode,  which  is 
an  electron  emitter,  has  some  of  the  same  characteristics  as 
the  anode,  but  is  more  complicated  to  represent 
analytically.  The  description  of  the  cathode  will  not  be 
attempted  here. 

Necessary  simplifications  limit  existing  probe  solutions 

to  special  cases  as  can  be  seen  by  sampling   the   literature 

regarding  plasma  probes.   Lara[10]  solves  a  sheath  problem  by 

matching  the  boundaries  of  "inner"  and  "outer"  solutions  and 

discusses   their   dimensionality.    Stahl  and  Su[11]  use  the 

same  approach  of  separating  the  sheath,  ambipolar  region  and 

free   stream.    Additionally,   they  prove  the  existence  of  a 

sheath  on  a  flat  probe.   Cohen[12]  criticizes  this   approach 

of   separate   regions   on   mathematical   grounds  saying  that 

quantities  are  forced  to  fit  and  some  derivatives  tend  to  be 

discontinuous.     McKee    and    Mitchner[13J    deal   with   a 

collisionless   sheath   (/^A.  )   but   include   ionization   and 

s 

recombination   in    the    ambipolar    region.    Bailey   and 

Touryan[14]  investigate  a  sheath  that  is  large  enough  to   be 

of   the   same   order  of  magnitude  as  the  boundary  layer,  and 

take  advantage  of  Blasius'  similitude  co-ordinates  to  reduce 

the  problem  to  one  dimension. 

Several  solutions  are  available  for  spherical  probes 
including  Kiel[ 15  J  ,  Barad  and  Cohen[16]  ,  Su  and  Lam[17]  , 
and  Cohen[12].  Kiel  ignored  diffusion  while  Barad  and  Cohen 
neglected  ionization  and  recombination.  Su  and  Lam,  and 
Cohen  were  the  first  to  use  a  systematic  analysis  of  probes 
in  collision-dominated  plasmas.  Lengyel[18],  by  dropping 
diffusion,  was  able  to  modify  the  elliptic  equations  in  such 
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a   way   as  to  solve  thera  by  the  method  of  characteristics  as 
one  might  do  with  hyperbolic  equations. 

Chung/  Talbot,  and  Touryan[19],  in  a  rather  extensive 
review  of  probe  work,  state  that  no  general  solution  is 
available  for  determining  charge  density  and  species 
temperature  for  probes  small  relative  to  boundary  layer 
thicknesses.  The  work  of  this  thesis  may  help  to  fill  that 
void  because  the  model  of  the  anode  is  that  of  an  array  of 
point  or  line  probes  immersed  in  the  non-convective  portion 
of  a  boundary  layer.  Results  of  this  work  may  also  be 
useful  to  other  fields  such  as  arc-discharges  and  lasers. 

C.   PHYSICAL  MODELING 

Both  theoretical  and  experimental  results  indicate  that 
in  the  MHD  environment,  where  pressures  are  near 
atmospheric,  the  sheath  lies  well  within  the  boundary  layer 
(about  10-5  m  for  the  sheath  thickness  or  three  orders  of 
magnitude  less  than  the  boundary  layer) .  Additionally,  the 
sheath  is  only  about  one  or  two  orders  of  magnitude  larger 
than  the  electron  mean-free-path  for  combustion  gases.  It 
becomes  apparent  then,  that  convection  plays  little  cr  no 
role  in  the  behavior  of  the  sheath  other  than  to  modify  the 
gross  temperatures  encountered  near  the  wall.  Consequently, 
the  boundary  layer  or  ohmic  problem  can  be  divorced  from  the 
sheath  and  treated  separately[ 20  ].  The  boundary  layer,  as 
well  as  the  wall  conditions,  determines  the  gas  temperature 
for  the  sheath  region.  Thus,  the  boundary  layer  voltage 
drop  can  be  added  to  the  sheath  drop  to  give  the  total  loss 
due  to  electrode  effects. 

As  mentioned  earlier,  analyses  of  boundary  layer  effects 
in  MHD  have  teen  developed,  but  they  tend  tc  be  complicated 
and   difficult   to  use.   In  order  to  have  a  complete  picture 
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of  electrode  drops  and  to  assist  in  interpreting  data,  a 
simplified  method  for  determining  the  voltage  drop  across  a 
thermal  boundary  layer  has  been  devised.  It  has  the 
capability  cf  rendering  quick  results  with  little  computer 
space  and  time.  It  does  not  pretend  to  be  a  complete 
analysis,  but  through  appropriate  assumptions  and 
simplifications  it  gives  the  drop  as  a  function  of 
temperature  only.  Appendix  A  gives  the  details  of  the 
method  along  with  its  limitations. 

An  investigation  of  the  nature  of  the  controlling 
eguations  for  the  sheath  shows  that  a  one-dimensional 
solution  does  not  exist  except  for  very  special  cases. 
This  investigation  resulted  in  the  publication  of  Sef.  21, 
and  a  summary  of  that  paper  may  be  found  in  Appendix  E.  It 
provides  a  basis  for  many  of  the  assumptions  used  in  this 
work.  The  following  conclusions  which  appear  in  Eef.  21 
helped  to  shape  the  model  chosen  for  this  work:  Current 
density  must  decrease  away  from  the  electrode  by  whatever 
means  possible.  The  mechanisms  of  geometry,  current 
constrictions  at  the  electrode,  or  ionization/recombination 
can  effect  this  decrease  singly  or  jointly.  A 
two-dimensional  flat  plate  in  the  absence  of 
ionization/recombination  would  be  at  least  partially 
one-dimensional  away  from  the  edges  and  is  therefore  not  a 
workable  model. 

Since  previous  solutions  have  been  obtained  with 
spherical  probes,  it  appears  that  a  geometrical  decrease  of 
current  density  is  a  proper  means  of  satisfying  the  nature 
of  the  equations.  Chen[22]  used  ionization/recombination  to 
extract  a  strictly  one-dimensional  solution  for  the  sheath. 
What  remains  then  is  to  study  the  features  that  result  from 
current  constrictions  at  discrete  periodic  sites  rather  than 
uniform  current  density  across  the  electrode  surface. 
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D.   OBJECTIVES  OF  THE  PRESENT  WORK 

In  this  work,  a  two-dimensional,  flat  plate, 
non-emitting  electrode  is  modeled.  Here,  the  current  has 
only  discrete  points  through  which  to  flow  on  the  electrode 
surface  (See  Fig.  1)  .  These  points  are  assumed  to  be 
periodic.  The  points  simulate  the  current  constriction  at 
the  electrode  which  may  be  there  because  of  temperature 
irregularities  or  surface  roughness.  A  collection  of  these 
points  or  nodes  might  represent  the  anode  spots  observed  by 
Biblarz[23]  and  Kimblin[24].  In  a  two-dimensional 
representation  the  point  is  actually  a  "wire"  of  unit  depth. 
The  spacing  cf  the  active  sites  must  be  greater  than  the 
sheath  thickness  in  order  to  avoid  one-dimensional  effects, 
and  small  enough  to  be  compatible  with  a  two-dimensional 
computational  field.  It  must  be  emphasized  here  that  this 
work  does  not  consider  the  phenomenon  of  arc  discharge, 
which  is  the  result  of  thermal  instaDilites[ 25 ].  Rather  it 
looks  at  the  current  constriction  required  to  satisfy 
continuity  of  charge,  Ohm's  law,  and  Poisson's  equation. 

The  primary  objective  of  this  work  is  to  present  a 
mathematical  model  of  a  two-dimensional  (Cartesian  geometry) 
electrode  in  a  quiescent  MHD  plasma.  The  controlling 
equations  are  solved  for  the  potential,  the  current,  and  the 
charge  density  distributions  in  tne  vicinity  of  the 
electrode.  The  additional  effects  of  a  magnetic  field  and 
Joule  heating  are  investigated.  It  appears  that  this  is  the 
first  investigation  of  the  magnetic  field  within  the  sheath. 
The  resulting  program  generates  sheath  and  ambipolar  regions 
in  a  self-consistent  way  using  the  same  set  of  equations 
throughout  the  field,  obviating  the  need  to  match  layers. 
The  size  of  the  sheath  and  the  voltage  drop  attributable  to 
its  shielding  effects  are  investigated. 
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E.   ORGANIZATION 

In  Section  II  the  physical  concepts  and  the  controlling 
and  boundary  equations  are  developed.  The  equations  are 
cast  into  a  form  for  use  in  the  computer  in  Section  III.  It 
includes  a  description  of  the  numerical  methods  which 
produced  successful  results,  as  well  as  a  short  chronology 
of  techniques  that  failed  to  work.  A  summary  of  the  results 
of  potential  and  charge  distributions,  and  current  and 
temperature  plots  for  all  the  cases  studied  is  included  in 
Section  IV.  Section  V  lists  the  conclusions  reached  from 
both  a  physical  and  numerical  point  of  view  and  gives 
examples  on  the  use  of  the  results. 

Appendix  A  shows  the  technigue  for  determining  tno 
boundary  layer  losses.  Appendix  3  is  a  condensation  of 
Ref.  21  proving  the  non-existence  of  a  one-dimensional 
solution  for  the  limit  of  the  non-reacting  flux  of  charges 
in  a  collision-dominated  plasma.  A  dimensional  analysis 
which  derives  a  set  of  independent  non-dimensional  variables 
for  the  controlling  equations  is  found  in  Appendix  C. 
Additionally,  Appendix  C  contains  a  fractional  analysis  of 
these  and  ether  equations  to  estimate  the  significance  of 
certain  physical  effects.  Finally,  the  computer  programs 
used  in  this  work  are  included  in  the  section  "Computer 
Programs",  with  explanations  as  to  their  use. 
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II.   ANALYSIS  OF  THE  SHEATH  AND  AMBIPOLAR  REGIONS 


A.   THE  SHEATH  ENVIRONMENT 

A  shaath  is  a  non-neutral  region  which  lies  adjacent  to 
the  electrode  or  insulator  surface.  The  strict  definition 
of  a  plasma  requires  charge  neutrality,  thus  the  sheath  is 
not  part  of  the  plasma  since  a  space  charge  exists.  For 
example,  the  anode  is  positive,  which  means  that  it  attracts 
electrons  and  repels  positive  ions.  As  the  electrons 
collect  at  the  anode  surface,  the  anode  potential  is 
partially  shielded  from  the  rest  of  the  plasma  by  the 
space-charge  potential  drop. 

Between  the  sheath  and  the  free  stream  lies  the 
ambipolar  region.  As  the  name  implies,  it  has  an  equal 
number  of  positive  and  negative  charges.  In  this  transition 
region  the  electrons  are  slowed  by  heavier  ions,  and  strong 
concentration  gradients  exist.  The  electrodes  themselves 
are  usually  metallic,  and  may  be  pins,  wires,  or  flat  plates 
separated  by  ceramic  insulator  wall  segments  (segmented 
electrodes).  The  surface  of  the  electrodes,  even  when  highly 
polished  have  roughnesses  of  the  order  of  the  sheath 
thickness,  and  in  addition  exchange  material  with  the 
flowing  plasma,  further  adding  to  the  irregularities.  These 
irregularities  increase  the  tendency  for  the  current  to 
constrict  by  providing  active  sites  for  the  current  to  flow 
along  minimum  energy  paths. 

The  size  of  the  sheath  in  MHD  generators  is  about 
10-s  m,  putting  it  well  within  the  boundary  layer  whicn  is 
about  10~2  m.  The  flow  is  therefore  essentially  stagnated 
in  the  sheath  region.  Though  the  ambipolar  region  extends 
further   into   the   boundary   layer,   a   fractional  analysis 
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(Appendix  C)  on  convection  shows  that  the  characteristic 
length  for  convection  in  most  flows  is  much  greater  than  the 
conduction  or  diffusion  lengths  within  the  sheath  and 
ambipolar  regions.  Consequently,  convection  may  be  entirely 
assigned  to  the  ohmic  study  of  electrode  losses  (see 
Appendix  A).  The  boundary  layer  lies  within  the  undisturbed 
plasma.  That  is,  it  lies  beyond  the  ambipolar  region,  there 
is  charge  neutrality,  and  electron/ion  concentrations  are  as 
would  be  predicted  by  the  Saha  equation[26]  for  plasmas  at 
local  equilibrium. 

In  the  model  presented  here,  the   following   assumptions 
are  made: 

1)  Steady  state, 

2)  Chemical  equilibrium  in  the  boundary  layer,  but  frozen 
flow  in  the  sheath, 

3)  No  induced  magnetic  field,  i.e.,  low  Magnetic  Reynolds 
number, 

4)  Negligible  ion  slip, 

5)  J   <<  J    , 

i     e 

6)  No  continuum  radiation  losses  in  the  energy  equation, 

7)  No  ion  emission,  and 

8)  Neutral  and  ion  particle  temperature  is  T  ,  unaffected 

o 

by  Joule  heating  and  uniform   within   the   domain   of   the 
sheath. 


B.   CONTROLLING  EQUATIONS 

1 •   Basic  Equations 

For  some  MHD  plasmas  the  mean-free-path  is  small 
enough  that  the  sheath  can  be  treated  as  a  continuum.  A  set 
of  collisional  equations  is  then  used  to  describe  the 
sheath,   ambipolar,   and   adjacent   free  stream  regions.   No 
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matching  between  these  regions  is  necessary  since  they  are 
generated  in  a  self-consistent  way  within  the  computational 
array.  The  equations  used  are  a  combination  of  the 
conservation  equations  and  Maxwell's  equations. 
Specifically,  they  are  expressions  of  conservation  of 
charge,  momentum  and  energy,  as  well  as  Gauss*  law. 

A  special  representation  of  the  electron  momentum 
eguation  is  embodied  in  a  generalized  Ohm's  law  and  may  be 
written  [Ref.  26,  p.  363] 

Js  "  -ys"seV*  -  -JFT  (JsxB)  *  DseVns  '  (1) 

where  the  +  sign  applies  to  the  electrons  and  -  is  for  ions. 

The  first  term  on  the  right-hand-side  represents  conduction, 

the  second  is  due  to  the  magnetic  JXB  current  and  the   third 

terra    represents   diffusion.    The   term   /3   is   the   Hall 

s 

parameter  and  gives  the   relative   effect   of   the   magnetic 

field   in   the   plasma.   The  equation  for  species  continuity 

can  be  written  as 

V-J   =  n  e  *2* 

s     s 

Appendix   C   shows   that  the   characteristic    length    for 

ionization/recombination  is   much  larger  than  the  she3th  or 

ambipolar  region  lengths.  Under  these  conditions  Eg,  2 
becomes 

V-J   =  0  P) 

s 

Poisson's  equation  comes  from  Gauss'  law  for  the  divergence 
of  an  electric  field  and  is  given  in  potential  form  as 

V2*  =  -(n.-ne)  e/e0  <"> 


Another  equation  that  will  be  Useful  later  and  which 
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applies    under    local    equilibrium    is   the    Einstein   relation. 


D   /u      =   kT  /e 
s'  Ms  s' 


(5) 


Electrons  absorb  energy  at  the  rate  of  J  »E  pen 

e 

unit  volume  from   the   electric   field.    Because   of   their 

relatively   small   mass,   the   electrons   are  inefficient  In 

transfering  energy  to  the  other  particles.   For  that   reason 

much    of    this   energy   goes   into   raising   the   electron 

temperature  above  that  of  the  heavier  particles.   Collisions 

with   these   particles,   both  elastic  and  inelastic,  tend  to 

limit  the  energy  retained  by  the  electrons.   In  a   quiescent 

plasma   an   expression   accounting   for   the  electron  energy 

balance  is  (Kef.  26,  p.  240-243) 

6  m 

J  -E  =  3n  k(T  -T)   Y   v    -^ 
e       e   e      L        es   m 


(6) 


2 •   Hon-Diroensional  Parameters 

The  non-dimensional  parameters  developed  in  Appendix 
C  are  repeated  here  for  use  in  this  section. 


n.  =  n./n 
l     i   o 


LV 


n  =  n  /n 
e    e  o 


x  =  x/L 


<j>  =  <j>/cj> 


o 


y  =  y/L 


(7) 


v;here : 


<b      =  kT  /e 


L  =  e  /(e  kTo)    and 


n   =  (kT  e  /e2) 3 
o       o  o' 


T    is   the  temperature  of  the  neutral  gas  which  is  constant 
o 

throughout  the  region  of  interest,  and  independent  of   Joule 

heating. 
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Each   of    the   three   characteristic      parameters      $    ,      L 

o 

and   n    are   functions  only  of  a  characteristic  energy  kT  . 
o  o 

ft   now  represents  a  non-dimensional  electron   charge,   or   a 
e 

measure   of  the  deviation  from  the  equilibrium  density  value 

in  the  free  stream.   n\   is   a   non-dimensional   ion   charge 

density,   and   within   the   sheath  n  *n  ,  whereas  within  the 

/\      A  A       A  A" 

ambipolar  region  n  =n   and  n  <n      and  in  the   free   stream 

e   i       e   Saha 
plasma  n  =/i  =/i     .   In  fact,  for  purposes  of  this  work,  the 

e   i   Saha 
edge  of  the  sheath  will  be  where   n    approaches   n    within 

e  i     a 

1  %,      and   the  edge  of  the  ambipolar  region  will  be  where  n 

A  e 

approaches  n     within  \%. 
Saha 


3»   Ppisson1 s  Equation 

Using  the  non-dimensionalizing  scheme   on   Poissou's 
equation,  Eg.  3  becomes 


/\  o  a     *■*  A 

V  d>  =  ne  -  n. 


(B) 


**  •   Energy  Equation 


E  = 


Applying  the  non-dimensional  parameters  of  Eq.  7   to 
-^74>  gives  the  following  result: 


E  =  - 


(kT0) 


—  V<J) 


(9) 


In   the   current   equation,  account  must  be  taken  of 

the  effect  of  varying  electron   temperature.    Since   J  >>J 

e    i 

the   current  is  essentially  an  electron  current.   Neglecting 

magnetic  field  effects  and  using   Einstein's   relation,   the 

current  is  written  as 
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kT   T 
o 
With    the   introduction   of   a   diinensionless   temperature, 

e  =  T  /T  ,    the  non-dimensionalized  current  equation  becomes 
6  °       Ue(M0)5  e04 

J^  =  — 5 °  (11) 


where 


j  =  -n  V<b  +  6Vn  l  ' 

J     e  y      e 


Equation  11  will  contain  the  term  (1+/32)  in  the  denominator 
of  the  coefficient  when  magnetic  field  effects  are  included. 
The  dot  product  of  Egs.  9  and  11  yields 

y  (kT  ) 7e  5 

J  «E  =  -  t"t D*Vd)  *'->' 

e  11 


Non-dimensionalizing  n   and  letting  d  =  ^V    0  Hie.  in   Eg.   6, 

e  s  es  s  rn 

and  introducing  Eq.  13  into  Eg.  6  results  in 

U     (kT    ) 3e    2 
<*  o  <~>  ^  -^ 

j.V<J>   =   n    (6-1)  (14) 


_       5  J       Y  e 

3ae 


As  the  electron  temperature  increases,  the  collision 
frequency  with  neutral  particles  increases  affecting  the 
terms  jJL        and  Q  .   Since  electron-neutral  collisions  are  by 

far  the  most  numerous  Q  can  be  approximated  by  CL=  V   0  — *L 

en  n  mn 

°r  /8kT  * 

a  =  n  Q   6  (m  /ra  )  "V 0  1/2  (15) 

n  en  n  e'  n   v  ra 

e 

In   the   range   of  interest  of  the  parameters,  the  collision 
cross-section  changes  little  with  temperature  so  Ct      depends 

1/2 

on  e 


30 


A    modification   to  Compton's   equation[ 27 ]      for      small 
values    of    E/P   gives    an    electron    mobility    of    the    form 

y     -  e9  ~1/2 (16) 

e  /8kT^ 

me  -\J7UT     Qennn 
e 

-1/2 

so  that  U      gees  as  6     .   Combining  Egs.  14 ,    15,  and  16  and 
*e 

using  the  perfect  gas  law  results  in  the  form  of  the   energy 

equation  used  in  this  work. 

,2  QV  (17) 


where 


YJ"V<J)  =  -n  (6  -6) 

_  0.13(kTQ)4  £q2  mn  (18) 

Y  =   2„  2   .   4 
p  Q   m  o  e 
r   en  e  n 

The  value  of  the  coefficient  "V  appears  to  dictate 
the  degree  of  Joule  heating  in  the  system.  This  can  best  be 
seen  by  expanding  Eg.  17  and  expressing  6  iraplicitely  in  the 
form 

(Y/6)VS-V$  +  1 


/\   /S    S\  /S  A 


(y/0)V(})-Vne/ne  +  1 


(19) 


For  physical  reasons  it  is  expected  that  6>1.  It 
can  be  seen  from  Eg.  19  that  for  very  small  values  of  ~y, 
6 — tl ,  which  is  the  constant  temperature  case.  Similarly, 
for  y  which  is  very  large,  an  upper  limit  for  6  is  found. 


•\  S\  /\  /\  /\ 


n    VcJ>-V(!>  (20) 

0    =    _1 for  y->oo 

V$-  Vn 
Y        e 

It   is   interesting  to  note  that  for  very  large  *y,  e  becomes 

independent  of  7,  which  means  that  it  is  independent  of   the 

chemistry  of  the  plasma. 


Now  that  y   has  been  described,  it  will  be  useful   to 
see   what   various   MHD   plasma   compositions  it  represents. 
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Assuming  atmospheric  pressure,  and  the  perfect  gas  law,   Sq. 

18  can  be  shown  to  be 

4 

AC  T   Z 

Y  -  10~46  x  -°  (21) 

O   Q 

nyen 

in  MKS  units  where  Z  is  the  atomic  weight  of  the  principal 
neutral  species,  or  a  mass  average  of  the  neutral 
constituents. 

For   combustion   type   gases   where   the    principal 

constituent  might  be  CO  ,  a  typical  value  of  7  can  be  found. 

2 

From  Bef.  26,  pages  139  and  148  for  a  temperature  of  1000°K, 

Q      =    1.5x10-****,   5        =      2x103,   y  is   found   to   be 
en  n 

98.3.   This   type   of   gas   has  y   values   which   are   low, 

characteristic     of     plasmas     with     high     collision 

cross-sections. 

In   many   inert  gas  plasmas  nitrogen  is  used  to  help 

the   plasma   reach   equilibrium   more   quickly.    The    same 

reference   gives   typical  values  of  Q    =  5x10_20m2  and  5  = 

es  n 

8.   Then  Y=    1.4x10s.   This  is  already  at   the   upper   limit 

described  by  Eq.  20.  Even  higher  values  are  obtained  with 
argon  and  cesium  seed  to  the  order  of  109  because  of  the 
Eamsauer  effect[Eef.  27,  p.  31]. 


5 •   Species  Eg  uat ions  w it h  Magnetic   Field   Effects   and 
N^gligJ-ble  Joule  Heating 

When  Joule  heating  is  small  and  collisions  between 
electrons  and  neutral  particles  are  sufficient  for  efficient 
energy  exchange,  electron  temperatures  will  not  differ 
significantly  from  the  neutrals.  Furthermore,  ion 
temperatures  will  be  essentially  the  same  as  the  neutrals 
regardless   of   Joule   heating   because  of  their  large  mass. 
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That  is,  the  system  of  equations  is  stable  even   when   T  >T 

e   o 

because   it  has  been  assumed  that  T  =T  =  a  constant  which  is 

o   i 

unaffected   ty   Joule   heating.    The   following   expression 

replaces  the  energy  equation  for  this  case: 


T   =  T.  =  T 
e     1     o 


(22) 


Now   Eg.   1   is   solved   for   J   and  J   in  a  two-dimensional 

x       y 

Cartesian  system  with  x  and  y  coordinates. 

^+    De 


Jx   -    (iVr1    [-Vsnse 


-    6    (-u   n   e  •—•  ±    D   e  -r — 
Msv    Ms   s      9y  s      3y 


Jy  =  (l+gV1  [-P£nse  ||  ±   D 


Sx 

D  e 
<- 

3n. 


)] 


se  8y 


3n. 


+  3  (-u  n  e 
s    s  s  dx 


s  dx 


(23) 


(24) 


Substituting   Eqs.  23   and  24  into  Eg.  3  and  using  Egs.  5  and 
22  to  help  simplify  results  in 

-nsV2*  -  V^Vns  ±  (kTo/e)V2ns  +  Ss<^  |±  -  !^  |±)  =  0 

(25) 


Eguation   25   represents   two   equations,   one  for  electrons 
(+sign)  and  one  for  ions  (-sign) .   Now  Q     =    0  because  of  the 

high   ion   inertia.    Hence   the   two  equations  representing 
species  continuity  and  Ohm's  law  become 


-neV   <j)   -   V<j>-Vn      +    (kT  /e)V   n      +   3  (- 


,/  f±>   -  0  (26, 


•n.V2<j>  -    V4>-Vn.    -    (kT  /eJV2!^   =   0 


(27) 


Non-dimensionalizing   and   substituting  Eg.  8  into  26  and  27 
results  in 
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ne(fi±-fie)  -  ^'Vne  +  9\  +  3*  =  0  (28) 


V*W  ~  ^?'^i  -  $% =  °  (29) 


where 


/S  /\  /N  ^ 


9n   3d)    9n   3d)  /->,-t\ 


9x   9y    9y   9x 


6 .     Species    Egua tions    with   Non-Constant   Electron 
Temper at  are 

For   the   case   of  significant  Joule  heating  Egs.  28 

and  29  are  net   descriptive   of   the   species   concentration 

since   T  *T    everywhere  in  the  field.   It  can  be  shown  that 
e   o 

for  this  case  the  species  equations  become 

9n  (n.-n  )  -  0Vn  -V(J>  +  i  n  V9-V<j)  +  6zVzn   +  ±  9V9-Vn 
e   i   e       e   r    2   e     Y        e2       e 

0/   1-   9 J  99    1-   9 $99     .9$  8ne 

+  PC"  2  ne   -  ~  +  2  ne   -  ~  +   G  -  ~  (31) 

z     3y  9x         9x  9y      9y   9x 

fl    3(1)    3"e         1    fi    99    8£e    ,    1     ft99    8% 

9x   9y  9x   9y  9y    9x 

n.  (n.-n   )    -    Vn.-Vcf-   -    V   n.    =   0 
lie  1      Y  1 

Notice   that   for  9  =  1  and  \Jq    =    0  these  equations  reduce  to 
28  and  29. 


C.   BOUNDARY  CONDITIONS 

Many  boundary  conditions  can  be  hypothesized  for  wall, 
electrode,  free  stream,  and  upstream/downstream  locations. 
Free  stream  conditions  are  actually  those  at  the  edge  of  the 
ambipolar  region.  They  are  taken  to  be  zero  potential  and 
zero  space-charge,  and  a  Saha-predicted  total   charge.    The 
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electrode  node  potential  is  a  fixed  value.  Othei  wall 
boundary  conditions  are  not  so  easily  determined  since  they 
are  frequently  vague  and  not  well  identified  in  the 
literature.  Also,  since  this  model  is  new,  the  boundary 
conditions  used  must  be  compatible  with  the  system  that  is 
being  represented. 

Successful   runs  have  been  generated  using  the  following 

conditions  with  T   =  T  : 

e     o 

1)  Free  stream  as  mentioned  above, 

2)  Periodic  electrode  nodes  allowing  periodic  conditions 
at  upstream  and  downstream  interfaces  (Fig.  1), 

3)  Line  electrode  nodes  allowing  current  constrictions 
(the  necessary  condition  for  two  dimensions  [21]).  The 
potential  is  fixed  and  ion  and  electron  concentrations  are 
zero  at  the  node,  as  dictated  by  a  catalytic  surface, 

4)  The  inactive  portion  of  the  electrode  wall  is  treated 
as  an  "insulated  wall".  Current  perpendicular  to  this 
wall  is  zero,  and  the  wall  eliminates  space  charge  and 
stagnates  charge  motion. 

In  reality  a  small  sheath  will  form  along  the  insulated 
wall  because  of  the  difference  between  electron  and  ion 
mobilites.  The  analysis  of  Appendix  C  shows  that  the 
thickness  of  the  sheath  goes  approximate] y  as  the  square 
root  of  the  potential.  Since  the  floating  potential  will 
generally  be  much  less  than  the  electrode  node  potential 
drop  it  is  expected  that  the  insulated  wall  sheath  will  be 
very  much  smaller  in  size  than  the  electrode  sheath. 
Consequently,  the  hypothetical  boundary  conditions  at  the 
wall  specified  above  are  representative. 

A  closer  look  at  the  electrode  node  boundary  conditions 
is  needed.  Blue  and  Ingold[28]  point  out  that  though  most 
authors  use  zero  charge  density  at   the   electrode   surface, 
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some   cases   require   the   modification   n  =  2J/(ec)  for  the 

s 

electrode    boundary    conditions.     Appendix     C     gives 

characteristic   values   for  J  and  c  of  7x10*  amps/ro2  and  10s 

m/sec,  respectively,  at  typical  temperatures.   This   results 

in   n  =  1019   particles/m3  which  is  of  the  same  order  as  the 
s 

charge  densities  which  will  be   encountered   later   in   this 

work.    The   difficulty   in  the  application  of  this  boundary 

condition  comes  from  the  fact  that  a  true  "point'1   electrode 

node   would   have   an   infinite  current  density.   As  will  be 

seen  later  the  numerical  procedure  uses  a   finite-difference 

scheme  and  the  electrode  is  not  a  true  "point".   The  current 

density  at  that  location  becomes  artificially  dependent  upon 

the   mesh-size.    To   avoid   this  problem,  either  the  charge 

density  must  be  set  to  zero,  or  some  new  non-zero   restraint 

must   be   found.    The   use   of  periodic  boundary  conditions 

seems   indicated   since   solutions   of   Laplace's    equation 

(Poisson^   equation   in  the  absence  of  a  space  charge)  tend 

to  be  periodic  in  two  dimensions.   As  will  be  seen,   several 

conditions  were  tried  and  the  results  are  discussed  in  later 

sections. 

Other  boundary  conditions,  such  as  a  catalytic  insulated 
wail  were  attempted  with  no  reasonable  results.  This  is 
discussed  further  in  the  next  section. 

Figure  2  is  a  representation  of  the  domain  in  the 
vicinity  of  the  electrode,  and  some  boundary  conditions  that 
are  typical  of  those  used  in  this  work.  The  domain  is 
represented  numerically  by  a  two-dimensional  array  of 
equally  spaced  points  which  are  operated  on  by  finite 
differencing  the  controlling  equations. 


Boundary  conditions  for  %r    n  ,  and  n   remain  the  same  as 

e       i 
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above  for  the  Joule  heating  case.  Temperature  boundary 
conditions  are  straightforward.  End  conditions  are  again 
periodic  and  free  stream  reguires  that  9=1.  Since  both 
the  electrode  node  and  walls  consist  of  dense  materials  with 
which  electrons  collide,  it  is  assumed  that  the  electrons 
give  up  their  surplus  energy  readily  to  these  surfaces. 
Therefore,  0=1  along  that  boundary. 
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FREE   STREAM    BOUNDARY 
<£>    =   0 

BC.:  ne  =  nsaha 

ft i         H saha 


OPEN    BOUNDARY 
B.C.:  PERIODIC 


ELECTRODE    NODE 
<£>  =  CONST. 

n,  =ne=o 


ELECTRODE  WALL 

■'e        ''saha 
B.C.'-   ni    =nsaha 


J 


=  0 


Figure   2 


Arrangement  of  electrodes  and  boundary  conditions 
for  computational  array. 


38 


III. 


METHOD  OF  SOLUTION 


A.  INTRODUCTION 

The      governing    eguations    (namely    Egs.    8,    28,    and   29)     are 

solved      numerically      for      the      parameters      $,      n    ,       and        n 

e  i 

throughout   the   domain   as  shown  in  Fig.  2.   Later  when  the 

electron  energy  equation  (Eg.   17)   is   included,   a   fourth 

parameter   6   is   included  in  the  field,  with  Eqs.  31  and  32 

representing  the  species  eguations.   The  method  of   solution 

will   now   be  discussed  together  with  the  application  of  the 

boundary  conditions.   The   system   of   partial,   non-linear, 

second-order   differential   equations   is  basically  elliptic 

and  constitutes   a   boundary   value   problem.    The   partial 

differential   equations   are   replaced   by  finite  difference 

equations;  the  technique  for  solution   is   discussed   below. 

The   non-linear  nature  of  the  controlling  equations  requires 

rather  sophisticated  computer  techniques,  and   the   boundary 

conditions  present  special  problems. 

There  are  a  number  of  possible  approaches  to  the 
numerical  analysis  of  the  equations  and  several  were  tried 
before  a  successful  approach  evolved.  In  order  to  better 
understand  the  nature  of  these  complicated  equations,  this 
Section  also  discusses  the  unsuccessful  approaches  and 
probable  reasons  for  failure. 

B.  NON-LINEARITY  CONSIDERATIONS 

The  twc  species  equations  contain  the  non-linear  terms 
which  present  special  programming  problems.  One-step 
techniques  for  solving  this  system  of  three  simultaneous, 
second  order,  non-linear,  partial  differential  equctions  are 
non-existent.   One   solution   procedure   in  such  cases  is  to 
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include  all  non-linear  terms  on  the  "right  hand  side,"   that 

is,   external  to  the  coefficient  matrix,  and  hope  that  these 

non-linear  terms  change  slowly  enough  with  each  iteration  to 

render  a  convergent  process.   This  procedure  is  known  as  trie 

Jacobi  method[29].   Looking  at  eguations  28  and  29   foe   the 

constant   termperature   case,   or   31   and   32  for  the  Joule 

heating  case,  it  can  be  seen  that  only   one   term   in   each, 

2         2 
either  V  n    or  V  n  ,   is   linear,   and   that   there  is  no 
e         i 

prescribed  "load"  on  the  right  hand  side.   As  will   be   seen 

in   the   expansion   of   these  equations  later,  this  leaves  a 

dozen  or  more  non-linear  terras   for   the   right   hand   side. 

Experience   shows   that   the   Jacobi   method   proves   to   be 

unstable. 

As  a  consequence  of  the  failure  of  the  Jacobi  method,  a 
quasi-Jacobi  method  is  used.  When  the  product  of  two 
variables  is  encountered,  one  variable  is  treated  as  a 
constant  coefficient  for  each  iteration.  This  means  that 
the  non-linear  terms  are  retained  in  the  coefficient  matrix. 
The  "constant"  coefficients  are  updated  after  every 
iteration,  thus  changing  the  coefficient  matrix. 
Successively  higher  voltages  were  calculated  in  increments 
allowing  small  changes  to  the  system. 

Had  this  quasi-Jacobi  method  not  succeeded,  the  next 
step  would  gave  been  to  apply  a  Newton-Eaphsen  technique^ 29  ] 
to  the  equations.  Because  of  the  increase  in  the  number  of 
terms  and  the  consequent  increase  in  computer  requirements, 
the  Newton-Raphsen  method  was  reserved  as  a  last  resort. 

C.   GRID  SIZE  AND  COMPUTATIONAL  SEQUENCE 

The  number  of  points  required  to  define  a  domain 
sufficiently  large  to  show  the  sheath  and  ambipolar  regions 
is  critical  (See  the  section   on   numerical   results) .    For 
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example,  if  a  two-dimensional  field  consisting  of  a  square 
array  of  31  x  31  points  is  desired  to  be  solved 
simultaneously,  then  it  can  be  shown  that  a  finite 
difference  scheme  will  require  in  excess  of  66  million  bytes 
of  storage  for  the  coefficient  matrix  alone.  Special 
storage  techniques  which  take  into  account  the  banded 
structure  of  the  coefficient  matrix  reduces  the  storage 
requirements  to  12  million  bytes.  Even  this  is  about  12 
times  the  core  capacity  of  the  IBM  360/67  at  the  Naval 
Postgraduate  School.  Obviously,  to  have  a  sufficiently 
large  array  and  to  keep  within  the  core  capacity  the 
solution  method  must  consider  portions  of  the  array  at  one 
time. 

This  problem  is  met  by  using  the  line  iterative  method, 
i.e.  by  solving  one  line  of  the  array  at  a  time,  sweeping 
back  and  forth  with  successive  field  iterations.  By  doing 
so,  the  size  of  the  coefficient  matrix  is  reduced  to  a  small 
fraction  of  the  full  matrix.  Because  the  equations  are 
elliptic,  and  therefore  each  point  affects  every  other  point 
in  the  field,  it  is  desirable  to  reduce  the  computer  time 
for  "information"  to  move  through  the  field.  This  is 
accomplished  by  rotating  the  field  90°  after  every  second 
field  iteration  such  that  the  line  sweeps  go  forward,  then 
backward,  then  up  and  then  down.  This  approach  takes 
considerably  less  execution  time  for  a  converged  solution 
than  the  usual  procedure  without  alternating  directions  and 
rotation . 

The  final  grid  size  chosen  for  this  line  iterative 
method  is  51  x  51.  This  decision  was  based  upon  the  time 
required  to  reach  a  converged  solution.  The  time  for  each 
run  was  of  the  order  of  four  or  five  hours.  The  coefficient 
matrix  requires  only  51  x  3  x  10  x  6  =  12,240  bytes  while 
the  total  space  needed  is  150,000  bytes  for  the  constant 
temperature  case,  and  170,000  bytes  for  Joule  heating.   This 
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storage    includes    program   space,   variable   storage   and 
solution  matrices  for  double  precision  arithmetic. 

To  further  increase  the  speed  of  convergence  a  method  of 

successive  over-relaxation  is  incorporated  similar  to   that 

outlined   in   Ref.   [30],   At  the  end  of  each  line  iteration 

the  parameters  are  subjected  to  the  equation 

A      =  (1  +  w)A     -  w  A  ,  ,  (33) 

corr  new      old 

where  A      is  the  corrected  value  of  the  parameter,  A     is 
corr  new 

the  most  recent  calculation,  A     is  the  previous  A     ,  and 

old  corr 

w  is  the  over-relaxation  parameter. 

D.   CONTROLLING  EQUATIONS 

The  finite  difference  representation  of  all  terms  was 
expressed  to  the  order  of  h2  when  practical,  where  h  is  the 
grid  spacing  of  the  square  51  x  51  array.  Ketter  and 
Prowell  (Ref.  29,  p.  226-228)  also  describe  the  procedure 
for  setting  up  finite  difference  equations.  Second 
derivatives  in  two-dimensions  are  given  by 

V^A.   .  =  (A...  .  +A.  ,  .  +A.  ,.n+A.   •  .-4A   .)/h2     (34) 
1,3      i+l,D    i-l'D    i/D+1    i/D-1     i/D 

First  derivatives  have  the  form 

tr-=(Ai+1,j-Vi,j>/<2h'  (35) 

1/D 
Where   a   first   derivative   is  required  at  a  boundary  it  is 
written  as  either  a  forward  or  backward  difference. 

^-=(-3Airj*4fli+1#j-Ai+2/j)/(2h) 


3xio 


Transforming    Eg.    8    into    finite   differences  gives 

b. ._    .  +  J.   .    .  +  $.    .._  +  $.    .  ,  -  4$.    .  +  h  n. 

^1+1,3       Ti-1,D       vi,3+l       yJ-,D-l  it!  \fj 


-  hi  n  =0 


U2 


This  is  the  simplest  of  the  three  controlling  equations. 

The  species  equations  are  considerably  more  complicated. 
The  finite  difference  form  of  Eq.  28  is  transformed  here  to 
illustrate  the  use  of  variables  as  "constant"  coefficients 
in  the  non-linear  terms. 


(*i,j+i  -  *i,j-i>  i^(no 


+    (n 


+ 


-  n  )    -   n 

e  „.    .,    ±  e 


+  n 


irj  +  1 


i/j-l 


i+l,j  i~l,j 

4n  )/h2 

e .     . 


i,j+l 


e . 


i,j-l 


/4h' 


A                              A                                    s\  s\                                              /\  /\  /\ 

n     (n     -  n.    )  +  ( <b .  .  ,  .  -  <h.  ,  .)  [n       -n 

e.  .   e.  .    i.  .  yi+l.n    yi-l.i  e.,n  e.  , 

if]    ifD     i/D  i+l, D  1-1/ 


+  $(n       -n      )  ]/4b/  -  (n 
ei,j+l   ei,j-l 


+  n      )/h    (38) 
ei+l,j    ei-l,j 


This  equation  corresponds  to  a  line  in  the  y  direction.  The 
bracketed  terms  {}  are  treated  as  constant  coefficients  but 
are  updated  with  each  field  iteration.  The  finite 
difference  version  of  Eg.  29  resembles  Eg.  38,  but  is 
considerably  simpler  since  /3.-0. 


When  the  electron  temperature  is  uniform   and   equal   to 

the   neutral   temperature,   no   energy   eguation   per   se  is 

reguired  since  the  condition  T   =  T   is  incorporated  in   the 

e    o 

two   species   equations.    However,   for  the  non-equilibrium 

case,  the  energy  equation   is   solved   separately   from   the 

other   tnree.    That   is,   Egs.   8,   31,   and   32  are  solved 

simultaneously  for  assumed  values  of  9.   With  each  iteration 

0  is  updated  by  solving  Eg.  17. 

E.   BOUNDARY  CONDITIONS 


U3 


Free   stream   boundary  conditions  are  the  simplest.   For 

each  case  the  free  stream  values  are  set  at  $  =  0,  n   =  n   = 

e     i 

n     ,  and  0=1. 
Saha 


Upstream/downstream  conditions  are  chosen  to  be 
periodic.  Ihey  require  that  not  only  the  values  of  the 
unknowns  be  the  same  at  periodic  points,  but  also  that  their 
derivatives  be  the  same.  For  example,  if  the  far  upstream 
station  is  labeled  " 1 "  and  the  far  downstream  station  is  "a" 
the  two  equations,  then 


A     =  A,  .     and  , 
n,D     1,3 


(39) 


A,   .  =  A-  .  -  A   .  +  A   ,   . 
1/3     2,3  n,3     n-1,3 


(40) 


A       A. 


would  apply  for  all  parameters  §,  n  ,  n  ,  and  9. 

e    i 


The   insulated  wall  is  hypothesized  to  be  neutral  and  to 
equilibrate   the   charges,   therefore,   the   conditions   are 


numerically   set   to  n   =  n 

e     i 


The  wall  condition  is 


Saha 

taken  from  the  restriction  that  perpendicular  current  is 
zero.  From  a  non-dimensional  form  of  Eg.  24  for  zero 
current  in  the  y-diirection  the  boundary  condition  at  the 
wall  is 


34>    3n 


-n 


3y   3y 


34> 
$n  —  + 
3x 


3n 


3x 


=  0 


(41) 


The   x-derivative   of   n    is   zero   along  the  wall  and  n 

e  e 


Saha 


So 


34>    3n 


saha 


3? 


34) 

T    Pnsaha  T"X 
3Y  3x 


(42) 


The  finite  difference  expression  for  Eg.  42  is 
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n   ,  6.  .  -  4n   ,  S.  .,.  -  3n   .   +  4n 
sahavi/j      sahaYi,j+l      saha     e.  .,, 

i/3  +  l 

=  -»saha*i,j+2  +  %  ,.,  +  *fisaha(?i+l, j  "  Vl, jJ 


(43) 


The  boundary  conditions  at  the  electrode  noda  consist  of 
three  conditions  $  =  <f> 

with  catalytic  electrodes. 


the  three  conditions  <J>  =  $  ,  n   =0,  and  n   =   0   consistent 

o    i  e 


F.   EVOLUTION  OF  THE  MODEL  AND  PHIOR  ATTEMPTS 

The  factors  affecting  the  choice  of  the  calculation 
procedure  and  computational  model  are  many.  Boundary 
conditions,  array  size,  electrode  placement,  choice  of 
coordinate  system,  and  the  form  of  the  controlling  equations 
are  but  a  few  of  the  things  that  had  to  be  considered  in 
modeling  the  problem.  The  final  successful  technique 
evolved  by  trial  and  error.  The  following  paragraphs  give  a 
brief  account  of  the  evolution  including  successes  and 
failures. 

At  the  outset  a  two-dimensional,  flat  plate,  continuous 
electrode  model  was  proposed.  Because  of  the  small 
thickness  of  the  sheath  relative  to  the  electrode  size,  most 
areas  of  the  sheath  would  not  experience  end  effects,  and 
would  be  effectively  one-dimensional.  Further  investigation 
led  to  the  proof  of  the  non-existence  of  a  one-dimensional 
solution  presented  in  Appendix  B. 

It  then  became  obvious  that  the  key  to  this  problem  was 
in  two-  or  three-dimensional  current  constriction  since  the 
one-dimensional  Cartesian  geometry  offers  no  natural 
geometric  means  to  decrease  current  density  away  from  trie 
electrode.  The  electrode  is  then  represented  by  a  series  of 
nodes  as  shown  in  Fig.  1.   An  examination  of  the   literature 
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showed  that  this  mode  of  current  constriction  was  known  and 
even  helps  explain  the  electrode  spots  which  may  be  present 
on  the  anode. 

At  first  the  model  consisted  of  both  the  anode  and 
cathode  on  opposing  sides  of  a  field  of  computational 
points,  but  it  was  realized  that  a  field  of  this  size  could 
not  possibly  show  all  the  details  of  the  sheath. 
Conversely,  a  field  capable  of  showing  the  required  details 
could  not  contain  both  the  electrodes.  So  the  final  model 
consisted  of  a  series  of  line  electrodes  (for  simplicity, 
the  anode)  along  an  "insulated"  wall  in  a  periodic  field. 
The  opposite  boundary  would  then  represent  a  free  stream 
conditon.  The  other  two  boundaries  are  upstream/downstream 
boundaries. 

Once   the   model   was   chosen,   the   job   of  solving  the 

equations  was  undertaken.   Initially,  the  three   controlling 

equations   (the  energy  equation  had  not  yet  been  considered) 

were  solved  separately  for  <$> ,    n  ,  and  n  ,   respectively,   in 

e        i 

iterative   fashion.    The   equations   were  Eqs.  3,  26  and  27 

with  a  Hall  parameter  of  zero.    The   calculations   diverged 

almost   immediately  because  in  the  vicinity  of  the  electrode 

many  derivatives  are  large   and   values   of   the   parameters 

changed  too  rapidly  from  the  initial  "guessed"  values. 

It  became  apparent  then  that  the  three  equations  would 
need  to  be  solved  simultaneously.  Section  III.C  explains 
some  of  the  problems  of  such  a  technique  in  terms  of  overall 
storage  requirements.  In  the  light  of  this,  the  line 
iterative  method  was  chosen. 

Armed  with  a  simple  set  of  variables  and  the 
line- iterative  method  a  solution  was  attempted.  Convergence 
was  achieved,  but  very  slowly,  since  it  took  many  iterations 
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for  "information"  to  cross  the  array.  So  the  array  was 
rotated  90°  at  intervals  causing  the  effective  sweep  to  be 
back  and  forth,  then  up  and  down. 

At  one  time  a  coordinate  transformation  was  effected 
putting  three  of  the  boundaries  at  infinity.  {The 
periodicity  cf  the  nodes  was  not  yet  considered.)  This  gave 
no  satisfactory  solution,  and  the  reason  is  believed  to  be 
the  following:  Although  Eg.  8  is  Poisson's  equation,  it 
reduces  to  Laplace's  equation  in  the  ambipolar  and  free 
stream  regions.  Since  there  is  no  solution  in 
two-dimensions  for  Laplace's  equation  satisfying  the 
requirement  that  <$>  =  constant  at  infinity,  th^n  no  solution 
to  this  problem  was  forthcoming. 

A   variable   mesh   size   was   also   attempted  with  some 

success,  but  since  it   was   determined   that   the  ambipolar 

region  was  of  the  same  order  of  magnitude  as  the  sheath,  it 
was  rejected  as  an  unnecessary  complication. 

Various  boundary  conditions  were  tried,  but  most  either 
did  not  work  or  gave  physically  meaningless  results.  The 
most  difficult  boundary  to  model  is  the  insulated  wall. 
From  the  final  boundary  conditions  chosen  the  system  of 
equations  describes  a  certain  physical  case. 

Once  the  present  technique  was  successful  for  various 
voltages,  it  proved  its  versatility  by  accepting  the 
addition  of  a  magnetic  field  and  the  energy  equation 
routinely.  While  not  being  the  ultimate  in  possible 
techniques  for  solving  this  system  of  equations,  the  present 
scheme  has  rendered  useful  results  and  promises  to  be 
further  productive  in  other  cases. 
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IV.   RESULTS 


A.   PROCEDURE  AND  CONVERGENCE  CRITERIA 

The  simplest  case,  which  included  no  magnetic  field 
effects  and  negligible  Joule  heating  began  with  the 
following  assumed  values: 

1.  zero  potential  everywhere  except  at  the  electrode, 

2.  selected  non-dimensional  potential  at  the  electrode, 

3.  zero    charge    density    everywhere    except    at   the 
boundaries,  and 

4.  Saha  charge  density  at  the  boundaries. 

Using  these  assumed  values,  iterations  for  improved  values 
proceeded. 

Periodically  convergence  was  checked  after  many  field 
iterations  by  comparing  the  left  hand  side  of  Eq.  0  with  the 
right  hand  side  at  various -locations  in  t^ie  field.  This 
convergence  criterion  may  appear  somewhat  arbitrary,  but  was 
chosen  when  it  was  found  that  Eg.  8  was  more  difficult  to 
satisfy  than  Egs.  28  and  29  and  therefore  represented 
stricter  convergence  criterion.  Additionally,  Eg.  8 
contains  all  three  unknowns  providing  a  check  of  a3 1 
parameters. 

The  error  was  calculated  from  the  finite  difference 
representation  of  Eg.   37. 

e  =  1     3-f3  +  1     i/P-1    yi+l>3    yi-l/3      1/3 

h  (n     -n.    )  (44) 

As  the  calculations  converged  the  error  decreased  with 
subseguent  iterations,  but  the  rate  of  convergence  slowed  as 
the   numoer   of  iterations  became  large.   Eventually  a  point 
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was  reached  where  further  iterations  were  not  cost  effective 
as  the  additional  improvement  in  the  solution  reguired 
excessive  computer  time.  A  second  convergence  criterion  was 
used  in  addition  to  the  the  first.  Certain  potential  values 
at  critical  points  in  the  array  were  plotted  witn  successive 
field  iterations.  If  these  values  approached  some 
assymptotic  value,  then  the  procedure  was  considered 
convergent. 

As  discussed  earlier,  a  successive  over-relaxation 
technigue  was  used  to  increase  the  rate  of  convergence.  It 
was  found  that  over-relaxation  parameters  exceeding  1.1 
caused  tho  field  to  diverge  for  all  but  very  low  voltages. 
Some  improvement  of  convergence  rate  was  noted  for  u  =  1.05 
so  this  value  was  used  for  most  runs.  When  attempting  to 
achieve  solutions  for  higher  voltages  or  other  conditions 
which  might  prove  divergent,  such  as  large  Hall  parameters 
or  Joule  heating  parameters,  a  value  of  w  =  0  was  used  in 
order  to  put  the  least  possible  "stress"  on  tha  system  of 
equations. 

B.   NEGLIGIBIE  JOULE  HEATING  AND  NO  MAGNETIC  FIELD  EFFECTS 

Convergence  for  this  set  of  runs  was  slow.  In  order  to 
achieve  an  error  of  less  than  20%  (|£|<-2)  over  1000 
iterations  were  required,  even  with  a  near  optimum 
over-relaxation  parameter  of  1.05.  This  equates  to  about 
four  hours  of  computer  on  the  IBM  360/67.  An  additional 
hour  (250  iterations)  reduced  the  error  to  16%.  Final 
results  shown  in  Figs.  3-18  were  completed  with  |£|<.16. 
Further  improvement  would  have  required  excessive  computer 
time. 

Figures  3,  6,    9,  and  12  are  plots  for  the  case  <t>   =  5.80. 
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This  value  represents  one  volt  at  T  =2000°K.   Figures  4,   7, 

o 

10,   and  13  are  for  $  =  17.41  and  Figs.  5,  8,  11,  and  14  are 

for  $  =  29.02,  representing  3  and  5  volts,  respectively,   at 

T   =  2000°K. 
o 


The  overall  size  of  the  array  for  each  of  these  runs  is 
500  times  the  characteristic  length  defined  in  Eg.  7.  At 
2000°K  this  means  that  the  dimensions  of  the  sgiiare  field 
shown  in  the  figures  are  5.25x10-5m  on  a  side.  The 
non-dimensional  charge  density  as  determined  by  the  boundary 
condition  is  10~3.  This  corresponds  to  8.5x1017 
particles/m3 ,  which  is  typical  of  an  equilibrium  charge 
density  for  seeded  plasmas  at  this  temperature. 

1  •   E5t£Ilii§.I  Distributions 

Figures  3-5  are  all  one-dimensional  profiles  of 
potential  and  charge  density  distributions  along  a  line 
perpendicular  to  the  electrode  wall.  The  aoscissa  shows 
distance  frcra  the  electrode  in  terms  of  the  characteristic 
length.  Figures  6-8  are  two-dimensional  contour  plots  of 
the  potentials  for  the  three  voltages.  The  electrode  is 
located  at  the  bottom  center  where  the  contours  concentrate. 
Except  for  the  absolute  magnitudes,  there  is  little  change 
noted  in  the  potential  plots  for  the  three  voltages.  As 
might  be  expected,  the  steepest  potential  slope  occurs  near 
the  electrode  where  the  greatest  electric  field  exists.  At 
the  free  boundary  the  slope  of  the  potential  is  nearly 
constant  indicating  a  constant  electric  field.  A  constant 
electric  field  is  a  property  of  two  opposing  infinite 
flat-plate  electrodes.  The  plasma  spreads  the  effects  of  a 
point  or  line  electrode  so  that  the  free  stream  behaves  as 
if  the  electrodes  were  infinite  plates. 
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2 •   Charged  Pa  r ticle  Distributions 

Figures  9-11  demonstrate  considerable  change  for  the 
electron  densities  with  increasing  voltage  as  do  Figs.  12-14 
for  ion  densities. 

An  examination  of  the  charge  density  plots  in   Figs. 

3,   4,  and  5  reveal  several  interesting  facts.   First,  since 

the  edge  cf  the  sheath  is   the   line   of   points   where   the 

charge   densities   are   egual,   the   figures   put  the  sheath 

length  at  35X,  907o,  and  100?,  of   the   total   length   of   the 

field.    In   the   following  table,  the  characteristic  sheath 

length  as  calculated  from  Eg.  C.7  for  T   =  2000°K   is   shown, 

o 

next   to   the   computer   generated   sheath   length   from  the 

extimate  cf  the  convergence  of  the  charge  density   lines   in 

the  figures: 

$         4>         characteristic    computer 
o         o 

sheath  length     generated 

-5  -5 

5.80       1.0  volts  .800x10    m        1.84x10    in 

-5  -5 

17.41      3.0  volts  1.38x10    m  4.72x10    ra 

-5  -5 

29.02      5.0  volts  1.79x10    ra        5.25x10    m 

These   offer  order-of-magnitude  comparisons  at  best,  but  the 

predictions  of  the  characteristic  sheath  length  with  Eg.  C.7 

are   useful   in   estimating   the   size   of  the  field  for  the 

computer  program. 

Notice  that  the  charge  density  slopes  on  the  right 
of  Figs.  4  and  5  are  significantly  greater  than  zero.  This 
means  that  the  field  was  not  sufficiently  large  to  encompass 
the  entire  ambipolar  region.  Since  the  boundary  conditions 
at  the  free  stream  are  imposed  as  constants  they  are 
"forced"  conditions  and  probably  result  in  error  for  these 
cases.  One  result  of  this  is  the  erratic  behavior  of  the 
charge  density  close  to  the   electrode.    As   will   be   seen 
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later,   this   behavior   has   little   effect   on   the  current 
distribution  as  long  as  it  occurs  in  a  confined  region. 

It  is  evident  that  the  observation  of  100%  of  the 
field  for  the  sheath  length  in  Fig.  5  is  probably  low  since 
the  convergence  of  the  two  charge  lines  occurs  artificially 
at  the  boundary.  To  prove  this,  and  the  other  observations 
made  above,  a  run  was  made  for  the  $  =  29.02  case  doubling 
the  field  length  to  1000  times  the  characteristic  length,  or 
10.5x10-5  m.  Figure  15  shows  the  dramatic  result.  The 
erratic  behavior  of  the  charge  density  lines  near  the 
electrode  is  gone  and  the  slope  is  essentially  zero  near  the 
free  boundary.  The  sheath  length  is  measured  at  aoout 
3.94x10-5  m,  2.3  times  the  characteristic  predicted  value. 
It  is  interesting  to  note  that  this  is  about  the  same  ratio 
as  the  treasured/predicted  ratio  of  the  $  -  5.80  case. 
Figures  16-18  show  the  two-dimensional  plots  for  this 
%   =    29.02  case. 

C.   NON-CONSTANT  ELECTRON  TEMPERATURE 

With  the  addition  of  the  energy  equation  the  effect  of 
Joule  heating  on  the  system  was  studied.  The  value  Y  =  100 
was  introduced  first  to  simulate  combustion  gases.  The 
results  are  shown  in  Figs.  19-22  for  the  cases  <f>  -  5.80  and 
29.02.  In  spite  of  the  earlier  discussion  showing  that  the 
size  of  the  field  is  inadequate  for  the  higher  potential 
case,  it  is  used  here  to  avoid  the  complication  of  changing 
the  electrode  spacing  when  comparing  different  potentials. 

The  nost  cbvious  change  caused  by  the  heating  of  tlie 
electrons  is  the  changed  slope  of  the  potential  indicating 
less  shielding  for  the  Joule  heating  case.  This  is  further 
shown  in  Figs.  23  and  24  which  represent  $  =  5.80  and 
Y=  103.     The    charge   density   profiles   are   relatively 
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unaffected  by  the  addition  of  high  temperature  electrons,  in 
fact  there  is  no  observable  change  in  the  sheath  length. 

Figure  25  is  the  electron  temperature  profile  (T  /T   vs. 

e   o 

distance   frcm  electrode)  using  the  same  cut  as  the  previous 

one-dimensional  plots.   This  plot  shows  that  the  majority  of 

the   Joule   heating  occurs  very  close  to  the  electrode  where 

the  largest  electric  fields  exist.   There   is   no   noticible 

electron   temperature   rise  in  the  arabipolar  and  free  stream 

regions.   This  graph  is  typical  of  the  temperature   profiles 

for   the  various  cases  and  differs  only  in  the  height  of  the 

peak  near  the  electrode.   This  peak  appears  to  be  a  function 

of  the  potential  as  well  as  the  Joule  heating  parameter. 

From  the  numerical  standpoint,  the  addition  of  Joule 
heating  offers  a  benefit  not  easily  foreseen.  Convergence 
of  the  system  with  the  addition  of  the  energy  equation  was 
accelerated  such  that  errors  of  less  than  2%  were  achieved 
in  less  than  250  iterations.  The  reason  for  this 
improvement  is  not  yet  clear  but  its  advantages  are  obvious. 
The  inclusion  of  Joule  heating  is  a  better  approximation  for 
most  generators,  and  provides  the  bonus  of  reducing  the  time 
of  numerical  convergence. 

D:   EFFECT  Of  EQUILIBRIUM  DENSITY 

An  interesting  case  is  that  of  increasing  the 
equilibrium  charge  density  to  much  higher  values  by  changing 
the  boundary  value  of  the  electron  and  ion  charge  densities. 
Presented  here  are  the  results  of  setting  ft  =  1  at  the 
boundary.  At  a  temperature  of  2000°K  this  corresponds  to  a 
charge  density  of  1021  particles/m3 .  Though  not  a  common 
case,  it  serves  as  an  indicator  of  some  of  the  numerical 
properties  of  the  method,  and  a  good  comparison  for  the 
sheath  length  prediction. 
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Figure  25.      Electron  Temperature  Profile,   $=5.80, y=100, 
P=0y   n=10_:5,   array  size=500L,   51x51  grid. 
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At   T   =  2000OK   and   for   $  =  5.80   and   n  =  1   at   the 
o 

boundary,  the   predicted   sheath   length   from   Eg.   C.7   is 

2.5x10-7   m.    At   this   same  temperature  the  characteristic 

length  as  defined  in  Eg.  7  is   L  =  1.0x10~7  m.    Figures   26 

and   27   show   the   case   where   the   field  length  is  25L  or 

2.62x10~6  m.   As  expected  the  sheath  length  is  about  10%     of 

the  total  field. 

The  dotted  portion  of  Fig.  27  is  drawn  to  agree  with 
previous  arguments  on  the  electrode  boundary  conditions. 
Actually,  some  of  these  higher  density  runs  were  made  using 
a  different  boundary  condition,  namely  one  that  attempted  to 
impose  a  consistent  total  current  through  the  field. 
Though  erroneous,  it  was  found  that  this  particular  boundary 
condition  has  little  effect  on  the  rest  of  the  field,  and 
that  conditions  were  actually  determined  by  the 
non-conducting  wall  surface.  The  "kinks"  seen  close  to  the 
electrode  in  several  of  the  one-dimensional  plots  appear  to 
be  a  resolution  feature  of  the  numerics  which  tends  to 
override  the  effects  of  the  charge  density  boundary 
condition  at  the  electrode. 


E.    EFFECT   OF   THE   NUMBER   OF   ELEMENTS   IN  THE  ARRAY  AND 
ELECTRODE  POSITIONING 

Earlier  trials  used  a  41x41  array  instead  of  the  51x51 
array  used  later.  In  general  the  smaller  arrays  did  not 
encompass  the  entire  ambipolar  region  as  can  be  seen  in 
Figs.  28-30  which  are  for  n  =  1.0,  J  =  11.6,  and  array 
size  =  20L,  25L,  and  35L,  respectively,  in  41x41  arrays. 
They  seem  to  lack  sufficient  room  to  satisfy  the  equations 
even  after  varying  the  mesh  size. 

Figures  28-30  also  show  the  changes  that  occur  when  the 
mesh   size   is   changed.    The   sheath   length  appears  to  be 
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consistent  with  the  scaling,  but  the  ambipolar  slope  does 
not  appear  to  change,  indicating  the  inadequacy  of 
mesh-sizing  to  control  the  size  of  the  field  when  the  number 
of  elements  is  too  low. 

For  a  qualitative  look  at  a  multiple  electrode  plot, 
Figs.  31-33  are  provided.  This  scheme  allows  a  current 
density  comparison  of  this  run  with  one  of  lower  potential 
and  the  necessary  symmetry  in  the  absence  of  a  magnetic 
field.   Current  profiles  will  be  considered  later. 

F.   MAGNETIC  FIELD  EFFECTS 

The  introduction  of  a  magnetic  field  induces  the  most 
visible,  although  not  unexpected,  effects  of  all  the 
phenomena  discussed  so  far.  As  might  be  anticipated,  the 
symmetry  of  the  field  is  lost  as  shown  in  Figs.  34-43. 
Notice  however,  that  while  the  potential  shows  significant 
distortion,  the  charge  densities  do  not,  at  least  for  the 
case  of  $  =  5.80. 

The  angle  of  incidence  of  the  potential  contours  with 
the  insulated  wall  is  approximately  the  complement  of  the 
Hall  angle  (arctangent  Q)  as  can  be  seen  in  Fig.  35  (45°) 
and  Fig  43  (21°)  .  This  is  easily  understood  by  solving  Eg. 
24   with   J  =0  and  considering  only  conduction  effects.   The 

y 


result  is 


(45) 


This  equation  says  that  at  points  where  J  =0  and  diffusion 

y 

is  small,  the  slope  of  the  potential  lines  is  ~p.  In  the 
vicinity  of  the  electrode  where  diffusion  plays  a  more 
significant  part  the  slope  begins  to  deviate. 
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As  with  the  runs  in  the  absence  of  a  magnetic  field  the 
effects  of  the  presence  of  a  point  electrode  are  dampened  as 
the  free  boundary  is  approached.  The  magnetic  field 
produces  little  or  no  distortion  to  symmetry  outside  the 
sheath. 


G.   CURRENT  EBOFILES 

For  each  of  the  previous  runs  a  complete  two-dimensional 
current  distribution  can  be  made  using  the  potential  and 
charge  density  distributions.  The  current  lines  make  an 
angle  with  the  potential  contour  lines  approximately  equal 
to  the  Hall  angle  and  would  theoretically  be  exact  except 
for  the  effects  of  diffusion.  Figures  44  and  45  show  the 
current  streamlines  in  the  field  for  $  =  5.80  and  (3  -  0  and 
1,  respectively.  These  were  sketched  from  numerical  data 
from  the  "Current  Program"  which  gives  the  current  in  vector 
form.  The  streamlines  entering  the  field  from  the  right  are 
from  the  active  site  that  lies  periodically  to  the  right  of 
the  one  in  the  field. 

To  determine  the  total  current  flowing  through  the 
system  it  is  necessary  to  look  some  distance  from  the 
electrode.  Theoretically,  for  this  model  the  density  of  the 
current  will  vary  from  some  given  value  at  the  free  stream 
to  infinity  at  the  electrode,  however,  as  explained  for  the 
case  of  charge  density  a  "smearing"  takes  place  due  to  the 
numerical  approximation. 

In  terms  of  the  amount  of  information  available,  Figs. 
46  through  48  are  the  most  important  in  this  work.  They  are 
the  current-voltage  diagrams  for  the  charge  densities 
ft  =  10-3  and  n  =  1.  The  dimensionless  current  density  was 
taken  one  or  two  stations  from  the  free  boundary  (in  order 
to  avoid  numerical  errors  at  the   boundary)   and   represents 
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Figure  kh.      Current  Density  Contour  Plot,  q)=5.80,y^O,p^O, 
n=10"3,  array  size=500L,  51x51  grid. 
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Figure  k5 .      Current  Density  Contour  Plot,  $=5.80,y  =0,p=l, 
n=10~3,  array  size=500L,  51x51  grid. 
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the  macroscopic  current  that  would  be  observed  in  a 
two-dimensional  model.  The  non-dimensionalization  of  the 
current  was  consistent  with  Eqs.  7  and  11  and  contain  the 
term  (I+/32).  The  conductivity  was  introduced  to  replace  the 
mobility,  and  constant  terms  were  collected  to  give  the 
simplified  coefficient  shown  in  the  figures.  An 
interpretation  of  these  results  as  well  as  examples  of  their 
use  is  given  in  Section  V  as  conclusions. 
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v-   CONCLUSIONS  AND  RECOMEND ATIONS 

A.   PHYSICAL  CONCLUSIONS 

A  successful  model  of  the  sheath  evolved  from  the 
assumptions  of  steady  state,  frozen  flow,  and  uniform 
temperature  for  neutral  and  ion  particles.  Qualitatively, 
the  results  of  the  analysis  used  in  this  work  offer  much 
insight  into  voltage  drops  attributable  to  electrodes  in 
contact  with  an  MHD  plasma.  Quantitatively,  these  results 
are  useful  if  it  is  remembered  that  the  model  is 
two-dimensional  and  does  not  give  the  additional  degree  of 
freedom  for  current  expansion  afforded  by  three  dimensions. 
A  summary  of  the  more  basic  conclusions  which  were  drawn 
from  the  results  is  presented  here.  They  will  be  explained 
individually  later  in  this  section. 

1.  The  sheath  can  be  self-generated  from  a  consistent  set 
of  equations  with  or  without  the  use  of  an  energy 
eguation. 

2.  The  electrode  node  boundary  condition  for  charged 
particle  density  has  little  effect  upon  the  field,  while 
the  insulated  wall  boundary  condition  has  a  profound 
effect. 

3.  Current  constrictions  are  necessary  at  the  electrode 
to  satisfy  the  system  of  equations. 

4.  The  resulting  current-voltage  diagram  has  a  curvature 
consistent  with  theoretical  predictions  for  space-charge 
in  a  one-dimensional  flow  of  current  and  experimental 
evidence. 

5.  The  current  density,  for  a  given  potential  drop, 
varies  inversely  as  the  node  spacing. 

6.  The  conductivity  of  the  plasma  within  the  sheath  is 
critical  in  determining  the  sheath  voltage  drop,  and 
methods   to   increase   this   conductivity  will  result  in  a 
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decrease  of  both  sheath  and  boundary  layer  drops. 

1 •   Sheath  Formation  and  Charge  Density  Profiles 

The  results  have  shown  that  a  sheath  develops 
consistent  with  theory.  The  thickness  of  the  sheath  varies 
with  the  potential  and  is  approximately  that  predicted  from 
theory  on  Debye  shielding.  In  each  of  the  resulting  plots 
presented  in  the  previous  section  the  size  of  the  ambipolar 
region  is  essentially  of  the  order  of  magnitude  of  the 
sheath  length.  No  quantitative  conclusions  have  been  drawn 
as  to  the  relative  size  of  the  ambipolar  region  other  than 
to  notice  that  higher  voltages  showed  the  tendency  of  the 
ambipolar  region  to  increase  in  length  relative  to  the 
sheath.  It  is  very  possible  that  at  higher  voltages  the 
ambipolar  region  may  grow  to  about  the  size  of  the  boundary 
layer  thickness.  But  since  the  potential  drop  in  this 
region  is  low  compared  to  that  in  the  sheath  it  is  not  as 
significant  in  this  analysis.  The  addition  of  Joule  heating 
and  a  magnetic  field  did  not  appear  to  affect  the  sheath 
length.  For  the  parameters  used  in  this  investigation,  the 
difference  between  unlike  charge  densities  (space  charge) 
was  about  10S5  of  the  equilibrium  density. 

Seme  of  the  resulting  one-dimensional  profiles  shown 
in  the  previous  section  exhibit  erratic  behavior  in  the 
first  few  stations  past  the  electrode  node.  This  is 
believed  to  oe  caused  by  numerical  convergence  problems  in 
the  procedure  arising  from  steep  derivatives  near  the  node 
and  the  extreme  non-linearity  of  the  equations.  This 
appears  to  have  little  or  no  effect  on  the  rest  of  the  field 
since,  as  can  be  shown,  the  current  distributions  derived 
from  these  runs  are  consistent  with  those  from  smoother 
runs.  Furthermore,  charge  densities  other  than  zero  were 
tested   as   a   boundary   condition   at  the  electrode  with  no 
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change  in  the  overall  results.  It  seems  then,  that  other 
conditions  mast  be  more  dominant  in  the  behavior  cf  the 
potential  field. 

One  very  significant  factor  appears  to  be  the  choice 

for   the   insulated   boundary   condition.    For   example,   a 

catalytic   wall   was   attempted   by   setting  n  =n  =0,  but  no 

i   e 

reasonable  results  were  obtained.   The  biggest  problem   with 

the  catalytic  condition  appeared  to  be  the  non-corapa facility 

of   the   wall   condition   with   the   free   stream    boundary 

conditions,   specifically   the   free   stream   slopes  did  not 

approach  zero  regardless  of  the  mesh  size  or   characteristic 

length.    Density   gradients  occurring  at   a  catalytic  wall 

occur  in  a  small  region  when  compared  with  the   sheath,   and 

violate   the   assumption  of  a  continuous  fluid.   Numerically 

this  means  that  the  charge  densities  must  go   from   zero   to 

equilibrium   values   within   one  grid  space.   Physically  the 

majority  of  the  field  should  not  "see"  this  wall  anyway.   So 

the   conditions   chosen,  namely  that  of  an  equilibrium  wall, 

are   consistent   with   the   rest   of   the   field   and   yield 

meaningful  results. 


2 •   Cur rent  Density  Distributions 

The  current-voltage  diagrams  in  Figs.  46  through  48 
show  the  behavior  of  the  system  for  various  conditions. 
Figure  46  shows  that  the  non-dimensional  current  density  is 
inversely  proportional  to  the  electrode  spacing  since  the 
current  density  for  the  $  =  29.02  case  with  the  array 
size=1000L  is  approximately  half  that  with  the  array 
size=500L.  This  is  shown  even  more  vividly  in  Fig.  47  where 
the  introduction  of  multiple  electrodes  with  double  the  step 
size  (keeping  the  electrode  spacing  constant)  resulted  in  a 
single  current  line.  The  significance  of  this  dependancy  is 
that   the  voltage  drop  for  a  given  current  density  becomes  a 
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function  of  the  electrode  node  density,  and  is  thus  a 
problem  which  has  inputs  from  surface  physics  as  well  as  HHD 
conditions.  It  must  be  remembered  that  for  the  square 
arrays  used  in  this  work  the  electrode  node  spacings  are  of 
the  same  crder  as  the  sheath  length.  The  minimum  spacing 
must  be  of  this  order  to  avoid  one-dimensional  effects  that 
would  result  from  further  "crowding"  of  the  active  sites. 
Spacing  the  nodes  much  greater  than  a  sheath  length  makes 
computation  more  difficult  since  a  much  larger  array  would 
be  needed  to  describe  tne  effects.  Thus,  this  work  models  a 
minimum  spacing  resulting  in  a  near  maximum  current  density 
for  a  given  potential  drop.  On  an  actual  electrode,  the 
spacing  wculd  be  random  and  dependent  on  empirical   factors. 

Figure  48  is  the  current-voltage  diagram  for  the 
ft  =  1  case  using  a  51x51  array.  It  was  noticed  that  for  a 
given  potential,  the  current  density  for  the  ft =  2  line  was 
approxiaately  twice  that  of  the  p  =  1  line.  Using  this 
result  a.  P  =  2  line  was  sketched  in  Fig.  45,  since  no 
converged  solutions  were  available  for  this  line  at 
n  =  10-3. 

The  introduction  of  the  magnetic  field  moves  the 
line  to  the  right  as  shown  in  Fig.  46,  and  appears  to  cause 
a  smaller  potential  drop  for  a  given  current.  But  recall 
that  the  current  is  already  multiplied  by  the  term  (1+/32) 
due  to  the  effective  reduction  in  conductivity,  and  this 
partially  offsets  the  reduction  in  the  potential  drop. 

The  introduction  of  the  energy  equation  shows  that 
Joule  heating  increases  the  current  for  a  given  voltage,  but 
only  to  a  slight  degree.  Ihis  is  because  the  region  that  is 
affected  by  the  Joule  heating  is  small  and  very  close  to  the 
electrode,  well  within  the  sheath. 

In    the    absence    of    a    magnetic    field    the 
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current- voltage  diagram  tends  to  be  concave  upward,  a  fact 
that  can  be  predicted  with  the  following  arguments:  Cobine 
(Ref.  27,  p.  129)  shows  that  the  high  pressure  space-charge 
current  should  go  as  V2/y3  in  one  dimension.  If  y  is  taken 
to  be  the  sheath  length  in  Eq.  C.7  then  it  can  be  shown  that 
V  goes  as  J2.  Figures  46  and  47  are  nearly  parabolic  in 
appearance,  and  have  curvatures  consistent  with  this 
analysis. 

The  zero  voltage  intercept  produces  an  interesting 
problem  since  it  has  different  signs  for  Pigs.  46  and  47. 
It  would  be  expected  from  Eg.  1  that  for  zero  potential 
difference  in  the  absence  of  a  magnetic  field,  the  electron 
current  density  would  be  positive  because  of  the  diffusion 
term  and  the  positive  charge  gradient  for  the  existing 
boundary  conditions.  However,  recall  that  the  validity  of 
the  catalytic  electrode  depends  upon  the  fact  that  it  is 
charged.  In  the  absence  of  a  charged  electrode,  the  node 
should  appear  the  same  as  the  insulated  wall,  and  no  density 
gradient  should  exist.  Therefore,  no  special  significance 
can  be  placed  on  this  intercept  without  further 
investigation  of  the  boundary  conditions. 

Reference  [31]  reports  on  the  beneficial  effect  that 

increasing  the  seed  fraction  has  on  lowering   the   electrode 

voltage   drops.   From  the  behavior  of  the  electron  collision 

cross-sections      with      temperature,      the      optimum 

conductivity[ 32  ]  is  reached  at  seed  fractions  higher  than  in 

the  core  because  of   the   lower   total   temperature   at   the 

electrodes.   Heretofore,  the  only  benefit  would  seem  to  be  a 

higher  boundary  layer  conductivity,  but  as   shown   in   Figs. 

46-48,   the  local  conductivity  also  affects  the  sheath  drop. 

Hence,  a  higher  value  of  G       is   expected   to   decrease   the 

o 

voltage   loss   associated  with  both  the  sheath  and  ambipolar 
regions. 
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An  example  would  be  appropriate  here  to  show  the 
orders  of  magnitude  of  the  quantities  being  dealt  with. 
Suppose  it  was  desired  to  find  the  potential  drop  across  the 
sheath  in  a  plasma  where  the  temperature  is  180G°K,  the 
current  density  is  0.2  amps-cm-2,  free-stream  conductivity 
is  0. 10  mhos/m,  and  the  charge  density  is  known  to  be  of  the 
order  of  1018  particles/m3 .  Equation  7  shows  the 
characteristic  charge  density  to  be  7x1020  particles/m3. 
Thus  Pig.  46  would  be  the  appropriate  current-voltage 
diagram  with  n=10-3.  The  non-dimensional  current  density  in 
the  absence  of  a  magnetic  field  is  1.50x10-5  which  gives  a 
dimensionless  potential  drop  of  32.0.  This  corresponds  to 
4.9  volts  at  the  given  temperature. 

3 •   Comparisons  with  Ex£er ^ment 

Argyropoulos    et    al[33]    cite    the    use   of   a 

sophistocat ed  computer  solution  for  the  boundary  layer  drop, 

comparing   the   results   with   an  experiment  on  the  AYC0-&PL 

channel.   Their  results  show  the  anode  drop  to  be  101   volts 

due   to   the   boundary   layer.    They  do  not  consider  sheath 

effects.     Using    these    data,    namely,    T     =  2000°K, 

wall 

J  =  3.95x10*  amps/m2,  /3  =  1.02  and  a  chemistry  consisting  of 

toluene  and  oxygen  with  cesium  seed,  the  conductivity  at  the 

electrode  is  1.14  mhos/m.   Figure  46  shows  a  non-dimensional 

current  density  of  4.307x10~5.   Using  the  [2  -    1   line   this 

gives  $  =  35.6  for  an  actual  drop  of  6.13  volts  at  the  given 

temperature. 

Another  comparison  is  from  the  work  of  Sonju  and 
Teno[34],  For  this  experiment  they  used  the  same  chemical 
constituents  as  the  previous  one.  The  conductivity  at  the 
wall  temperature  of  1800<>K  is  0.30  mhos/m  and  the  other 
parameters  are  3  =  1.6x10*  araps/m2  and  ft  =  1.8.  Figure  16 
gives  j  =  1.7x10-*,  $  =  121.3,  and  a  sheath  voltage  drop   of 
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18.8  volts. 

Caution  must  be  exercised  when  making  these 
comparisons  since  several  factors  have  not  been  considered 
yet.  First,  tendencies  towards  non-equilibrium  would 
increase  the  conductivity  near  the  electrode  resulting  in  a 
smaller  potential  drop.  Secondly,  nothing  has  been  said 
about  the  node  density  of  the  experiments  as  compared  to  the 
assumptions  made  in  the  numerical  solutions,  namely  that  the 
distance  between  nodes  is  of  the  order  of  the  sheath  length. 
Should  the  actual  node  density  differ  from  this  it  will 
affect  the  current- voltage  characteristics.  Thirdly,  the 
use  of  a  two-dimensional  model  rather  than  a 
three-dimensional  one  may  tend  to  predict  different  voltage 
drops.   The  following  example  illustrates  these  limitations. 

An  experiment  by  Kessler  and  £ustis[6]  was  conducted 


using   ethanol   in   oxygen   with   1 


KOH 


Nitrogen 


was 


introduced   for   cooling   such  that  the  N  /0   ratio  was  0.5. 

2   2 

For  this  run  T     =  1685<>K,  T     =  2700<>K,  J=.  75x1 0  ^amps/m? , 
wall  core 

P =    1.5,   and   the   electrode   conductivity   is  0.56  mhos/m. 

Figure  46   gives   j  =  3.697x10~«,  I    =  490.2,   and   $  =  71.1 

volts.    This   is   a   higher  drop  than  the  measured  45  volts 

that  should  account  for  both  the  sheath  and   boundary   layer 

contributions. 

The  primary  difference  between  the  above  examples 
appears  to  be  the  conductivity  at  the  electrodes.  Since  the 
equilibrium  conductivity  has  such  a  strong  dependence  on 
temperature,  the  electrode  temperature  is  a  critical 
parameter  in  controlling  the  sheath  losses.  A  very  accurate 
comparison  can  only  be  made  using  a  plot  generated  from  the 
correct  charge  density  and  on  precise  knowledge  of  anode 
tern  perat ures. 
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It   is   common   practice   to   predict    the    losses 
according  to  the  empirical  formula 

V  =  A  +  B6      '  u\S) 

where  S  is  the  boundary  layer  thickness  and  A  and  B  are 
constants  to  be  determined  for  a  particular  experiment. 
Such  a  scheme  might  be  useful  here  if  A  is  the  sheath  drop 
and  B  is  the  boundary  layer  drop.  Using  the  comparisons 
above  and  from  Appendix  A  for  the  boundary  layer  drop,  the 
formula  wculd  look  like 

Sonju  and  Teno      V  =  18.8  +  1190(5 

Argyropoulos  et  al   V  =  6.1  +  11100 
where  $  i£  measured  in  meters  and  V  is  measured  in  volts. 

B.   NUMERICAL  CONCLUSIONS 

The  numerical  technique  is  discussed  in  detail  in 
Section  III.  Generally  it  consisted  of  a  finite  difference 
scheme  for  the  solution  of  a  system  of  three  (later  four) 
coupled,  non-linear,  partial,  second-order  differential 
equations.  The  line  iterative  method  of  solution  was  used. 
This  method  allowed  a  smaller  storage  requirement,  and  a 
choice  of  a  51x51  array  balanced  the  computer  time  to 
converge  with  a  mesh  size  sufficient  to  show  the  detail  of  a 
two-dimensional  sheath.  Convergence  required  four  to  six 
hours  of  computer  time  on  the  IBM  360/67  computer  for  each 
potential  and/or  Kail  parameter  chosen.  The  addition  of 
Joule  heating  reduced  the  convergence  time  to  one  hour 
because  of  the  damping  effect  of  the  temperature  term  near 
the  electrode. 

The  use  of  successive  over-relaxation  to  reduce  the 
convergence  time  met  with  only  marginal  success.  Because  of 
numerical  instabilites,  an  over-relaxation  parameter  greater 
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than  1.1  usually  caused  the  process  to  diverge. 
Consequently,  a  value  of  1.05  was  used  resulting  in  less 
than  a  10%  improvement  in  time. 

In  its  present  form,  the  computational  procedure 
outlined  in  this  work  is  limited  as  to  the  size  of  the 
parameters  describing  the  system.  For  example,  divergence 
resulted  for  non-dimensional  potentials  exceeding  35,  Hall 
parameters  greater  than  2,  and  Joule  heating  parameters 
greater  than  10*.  It  can  be  observed  from  fig.  28  that  the 
Hall  parameter  weights  the  relative"  size  of  the  axial 
derivatives  compared  to  the  cross-channel  derivatives, 
giving  rise  to  numerical  instabilities  for  excessively  large 
values  of  .  The  maximum  potential  attainable  seeias  to  be 
limited  by  the  size  of  the  array  since  a  value  of  $  =  35  was 
obtained  for  a  51x51  array  and  only  29  was  achieved  for  a 
41x41  array. 

C.   RECOMMENDATIONS  FOR  FURTHER  STUDY 

Since  the  principal  limitations  to  this  work  arises 
because  cf  the  large  storage  and  time  required  for 
calculations,  much  more  could  be  accomplished  with  improved 
numerical  techniques.  A  search  by  the  author  for  a  packaged 
subroutine  capaole  of  solving  the  computat ional  matrix  more 
efficiently  than  the  one  used  here  proved  futile.  Improved 
techniques  are  needed  which  will  accomodate  larger  fields  of 
computation. 

The  validity  of  these  results  is  limited  by  the  fact 
that  the  current  density  decreases  in  two  dimensions  rather 
than  three.  A  three-dimensional  scheme,  admittedly  a  major 
task  to  develop,  could  render  more  realistic  information  on 
the  effects  of  the  sheath  and  on  other  electrode  phenomena. 
One  possible  scheme  is  to  model  a   three-dimensional   system 
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by  using  an  axi-symmetric  cylindrical  geometry  for  the  point 
electrodes. 

Additional  physical  effects  involving  changes  in 
boundary  conditions  are  possible.  For  example,  boundary 
conditions  could  change  the  anode  to  a  cathode  by  changing 
the  sign  cf  the  charge  and  making  the  electrode  an  emitter 
of  electrons. [Ref.  27,  p.  106]  The  addition  of  convection 
into  the  system  of  equations  (should  it  be  deemed 
appropriate)  is  expected  to  accelerate  convergence  because 
it  would  ease  the  requirements  of  current  density  decreasing 
away  from  the  electrode. 

Other  effects  such  as  neutral  and  ion  heating  may  be 
incorporated  into  the  scheme  to  investigate  thermal 
instabilities  known  to  exist  for  such  a  case.  Conduction  of 
heat  in  the  solid  walls,  ion  slip,  and  frozen  flow  may  be 
added.  It  is  anticipated  that  frozen  flow  would  decrease  the 
importance  of  the  boundary  layer  drop  and  increase  the 
importance  of  the  sheath  drop. 

With  each  additional  effect  the  model  rises  in 
complexity  and  the  numerical  scheme  becomes  more  unwieldy; 
therefore,  only  the  principal  effects  should  be  reflected  in 
the  description  of  the  sheath. 
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APPENDIX  A 

A  SIMPLIFIED  TECHNIQUE  FOR  DETESWING  THE  BOUNDARY  lAYiR 
VOLTAGE  DROP  IN  MHD  GENERATORS 

As  a  complement  to  the  two-dimensional  sheath  problem, 
and  in  order  to  assist  in  understanding  the  overall  loss 
problem,  a  one-dimensional  boundary  layer  analysis  was 
investigated.  The  result  is  a  method  for  determining  th_e 
voltage  drop  across  the  thermal  boundary  layer  which  is 
simpler  than  most  methods  in  existence  today. 

The  technique  balances  sophistocation  with  ease  of 
application.  Because  of  its  lack  of  rigor,  many  simplifying 
assumptions  have  been  made  which  render  the  results 
inaccurate  for  certain  cases.  These  assumptions  are 
detailed  in  the  following  paragraphs,  and  experimental 
comparisons  are  shown. 

1 .  Geometry 

Figure  A1  is  a  schematic  of  the  voltage  drops  across  a 
typical  channel.  the  geometry  assumes  opposing  electrodes, 
but  does  not  preclude  the  possibility  of  diagonal 
construction.  This  analysis  looks  at  only  the  ohmic  portion 
of  the  voltage  profile  and  assumes  that  the  sheath  losses 
can  and  must  be  determined  separately.  The  boundary  layer 
may  be  fully  developed  and  uses  the  wall  conditions  as  one 
of  the  boundaries  for  ohmic  analysis. 

2.  Boundary  Lay_er  Drop 

Reconstructing   and  modifying  a  formulation  by  Rcsa[35], 
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the  average  resistivity  ((7-i)   across   the   channel   can   be 
written  as 


D 


(a"1) 


o 


1  A 


r-        6 


_1_  1 

°c   D 


/   <-7T>  d* 


D-6 
/ 


D 


•f 


/ 


u0 


o, 


(—)  dy 

a 


D-6 


ac 
{•£)    dy 


(A.1) 


where  6  is  the  thermal  boundary  layer  thickness  and  D  is  the 
channel  width  between  opposing  electrodes.  Assuming 
constant  conductivity  in  the  free  stream  region,  symmetric 
boundary  layers  at  opposing  walls,  and   non-dimensionalizing 


the  integration  variable,  the  equation  becomes 

1 


(a'1)    = 


C  L. 


1  +  2 


D 


0 


c  ,,y.    „  6 

0     0        D 


(A.  2) 


Eosa  [Eef.  32,  p.  75]  assumed  a  plasma  where  the  gas 
properties  vary  only  across  the  channel,  and  averaged  Ohm's 
law.   The  resulting  equations  for  current  are 


-  §  (Ey  +  U  *  B)  "  I  Jx 


(A.  3) 


J  =  a  E  +  BJ 
x      x   w   y 


(A.  U) 


where  G-   (7  (-:£--)  ~p2   and   all   terras   have   been   averaged 

across   the   channel   (y-direction) .   The  x-direction  is  the 

direction   of   flow.     Eliminating   J  ,   assuming  Q      is   a 

x 

constant  across  the  channel,  solve  for  E   to  get 

y 


E,,  =  UB  +  -Z.    (G  +  B  )  =  3Ev 
y         o  * 


(A.  5) 


The  actual  output  voltage  is 
D 


V 


out 


=  f    E  dy  =  E  D  =  UBD  +  -^  (G+32)  +  3E  D 
j  v      v  o  x 


(A.  6) 
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The   intercept   of  the  voltage  profile  with  the  channel  wall 
(see  figure  A1)  is  given  by 

D  =  UoBD  +  jX  (1  +  B2)D  +  BEXD  (1<7J 

y v^ 


Taking  the  difference  and  substituting   for   G   the   voltage 

loss  is 

___       j  D 

Vloss  =  (U0"U)BD  "  (1  +  B  )[ac(a  1)  -  1]  -£-         (A. 8) 

c 

This  loss  is  Ohmic  and  is  attributed  to  the  t*o  boundary 
layer  temperature  drops.  For  a  turbulent  boundary  layer  the 
first  tern  will  be  small  compared  to  the  second 


J  D 


Vbl  =  2d1-  <1  +  e2>  [^(a"1)    1]  (A. 9) 


c 


3 .   Conductivity  E x.p_r  j?§_sio n 

An  expression  will  now  be  derived  for  7=?  to  be  used   in 

Eg.   A. 9.   Hurwitz,  Kilb,  and  Sutton[36]  have  shown  that  for 

cases  where  the  current  is  entirely   cross-channel  (Faraday 

mode)  ,   the   conductivity   can  indeed  be  treated  as  a  scalar 
since   for   roost   such   channels   the   magnetic    field   is 

approximately   uniform  between  the  electrodes.   For  the  case 
of  scalar  conductivity  we  can  write[ 26  ],[ 37  ] 

0   =  nee2/[mec-enn(  J  xkQk)  ]  (A.  10) 

k 

The  subscript  k  denotes  the  various  species  present   in   the 

plasma,  and  x  =  n  /n  .   For  non-reacting  plasmas  the  sum  of 
k     k   n 

the   products   x  Q    will   be   relatively    independent    of 

k  k 

temperature.    We   will  assume  this  to  be  approximately  true 

for  reacting  plasmas.   Taking  the  ratio  of  the   conductivity 

at  a  point  to  that  at  free  stream,  Eg.  A. 10  leads  to 
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n   en 
eo  e  n 


n  c   n  . 
e  eo  nO 


(A. 11) 


Assuming  a  locally  Maxwellian  velocity  distribution, 
guasi-eguilibrium  (allowing  the  use  of  the  Saha  relation) , 
and  constant  pressure  across  the  boundary  layer,  Eg.  A. 11 
becomes 

a 


c  _  a-7/4   a(l/0-l) 


a 


(A. 12) 


where  CL and  6  =  T/T  . 

2kT„  o 


4.   Constant  Specific  Heat 


At  this  point  a  description  of  the  turbulent  boundary 
layer  is  needed.  For  the  purpose  of  this  analysis  we  will 
use  the  1/7th-power  law  and  Reynold's  anaiogy[ 38  ]  to 
describe  the  temperature  field.  An  extensive  and  detailed 
study  of  the  heat  transfer  problem  related  to  these  high 
temperatures  was  conducted  by  Brira[39],  but  in  the  interest 
of  simplicity  the  Prandtl  Number  is  assumed  to  be  one.    For 


constant  C  , 
P 


6 


T  -  T 


w, 


T  -  T  ' 
o   w 


=  ( 


1  - 


J±)  7 


•V 


(A. 13) 


so 


.(?)  = 


7U  ~  <Jw> 


7 


de, 


(A.  14) 


Now   combining   Egs.   A. 2   and   A. 12   and  changing  v?riables 
according  to  Eg.  A. 14,  the  result  is 

1  +  2i   /  ,-V«e°'V8-l)7'e"V6  d9 


c 


D 


w 


(1-0  ) 
w 


-  2  D 


(A. 15) 


Putting  Eg.  A. 15  into  A. 9  for  the  average  resistivity  gives 
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Vbl   =  V1  +  3    )6/°c 


frr  J 


(i  -  e  ) 

*-  w  w 


,-7/4    «(i/e-i)(         )6  d6 ._ 

W 


(A.  16) 


5 •      l2mconst an t    Specific   Heat 


For    a   plasma    with    a   non-constant    C      the    1/7th-power      law 

P 


is 


h  -h 

1  =    (  w    \  7 

6         vh     -h   ; 
o        w 


A  change  of  variables  would  then  give 


(A. 17) 


Vbl  =  Jy<1  +  3  )6/(Jc 


7C  T 
P  o 

-  (h  -h  ) 
o   w 
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/  e 


-7/4  ea(l/8-l)        6 

w 


1 


w 


(A.  18} 


Obviously  an  enthalpy  table  or  equation  is  needed  to  solve 
Eg.  A.  18.  However,  Mollier  diagrams  are  obtainable,  and 
once  the  entnalpy  relation  is  known,  the  method  shown  in 
this  analysis  retains  its  simplicity. 

6-   Weighting  Filiations  and  Voltage  Drop_  Parameters 

Figure  A2  is  a  graph  of  the  conductivity  from  Eg.  A. 11 
for  Toluene  and  pure  oxygen  with  2%  cesium  seed  at  a  core 
temperature  of  2600°K.  Plotted  also  is  the  conductivity 
from  a  more  sophistocated  program  provided  by  AVCO  Everett 
which  also  assumes  quasi-equiiibrium,  but  does  not  make  the 
other  simplifications  done  above.  The  agreement  at  2600°K 
results  from  the  "forced"  condition  in  the  simplified 
program. 

An   examination   of   Eqs.   A. 16   and  A. 18  shows  that  the 


99 


c 

0) 
O 

01 

ex 

§ 

4-1 


•H 

a 

01 

bO     • 

o  o 
o 

II 

3 


0) 

(-1 

4-J 

a) 


o 

H 


0) 
CO 


o  • 

•H  CO 

^  a; 

•u  u 

e  -5 


ex 
w 
o 

e 

CO 


So 

l«     4-1 
J-l    X) 


•H 

> 
•H 

o 
X) 
o 


o  u 


AllAllOnQNOD 


o 

< 

o 

CO 

0) 

~"~ 

M 

3 

M 

•H 

tn 

100 


effect  cf  changing  the  integration  variable  was   to   average 

the   resistivity   across  the  boundary  layer  in  "temperature 

space"  instead  of  "geometry  space".  In  each  case  this  change 

generated    a    weighting    function    which   weights    the 

contribution  of  the  resistivity  to  the  integral  according  to 

the   temperature.    Other   weighting   functions  derived  here 

are: 

7(9-ew)  6 
for  constant  C  :        H   =  

P  H-ew)7 

(h-hw)  & 

for  variable  C  :         H   =  7C  T 

p  h      p  o(h0-hw)  7 

Other    descriptions   are   possible,   and   often   desirable. 

Useful   descriptions   for   laminar    boundary   layers    are 

available  in  Cramer  and  Pai[40]  and  Kerrebrockf 4 1 ]. 

Figure  A3  shows  the  effect  that  the  weighting  function 
has  on  resistivity.  The  resistivity  (reciprocal  of 
conductivity)  is  plotted  first  from  the  AYCO  program  and  the 
approximate  value  from  Eq.  A. 12.  Next,  the  reacting  plasma 
weighting  function  times  the  differential  Ae  is  plotted. 
Lastly,  the  product  of  the  first  two  graphs  shows  that  the 
effect  is  to  reduce  the  error  considerably,  by  weighting 
more  heavily  those  values  near  the  core.  In  fact  the 
average  resistivity  from  Fig.  A3  (area  under  the  right  hand 
graph)  is  0.0254  for  the  AVCO  program,  and  0.0244  for  the 
simplified  program.  Obviously,  the  effect  of  the  boundary 
layer  shape  is  to  smooth  out  the  errors  of  the  approximate 
resistivity  calculations. 

From  Egs.  A. 16  and  A. 18  a  non-dimensional  boundary  layer 
voltage  parameter  can  be  identified. 

♦bl  "  Vblac/Uy(l  +  62)«]  <A-19> 

so  that 
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*bl=  /' 


0 


-7/4   a(l/e-l) 


w(9)  de  -  i 


(A. 20) 


w 


where  W{6)  is  the   weighting   function   appropriate   to   the 

boundary   layer   description.   Figures  A4a  and  A^b  are  plots 

of  <J>    vs.  wall   temperature   and   core   temperature   for   a 
bl 

constant   C    gas   and   Cesium   seed.   The   former   uses  the 
P 

1/7th-power  law  while  the  latter  uses  a  1/8th-power  law.   As 

might   be   expected   the   1/8th-power  law  predicts  less  of  a 

voltage  loss  because  the  region  of  lower  temperature   plasma 

is   thinner.   Figures  A5a  and  A5b  are  similar  graphs  fcr  the 

combustion  products  of  Toluene  in  pure  oxygen   and  2%    cesium 

seed  added. 

7«   Experj. mental  Comparisons, 

As  with  all  experimental  predictions  the  accuracy  of  the 
results  is  dependent  upon  the  accuracy  of  the  input  data  as 
well  as  the  ability  of  the  model  to  represent  the  actual 
situation.  The  following  examples  illustrate  the 
convenience  of  this  method,  but  at  the  same  time  show  the 
importance  of  an  adeguate  model. 

Sonju  and  Teno[34]  report  the  results  of  an  experimental 
run  on  the  Viking  I  generator.  Fig.  A6  is  a  reproduction  of 
the  voltage  profile  for  one  run.  Because  of  burner 
innef f iciencies  they  found  the  conductivity  of  the 
Toluene/oxyqen/cesiura  plasma  to  be  between  12  and  16  mhos/m, 
about  half  the  33  mhos/m  that  would  have  been  predicted   for 


complete    combustion[  42  ].    Using  {3   =1.8,   J 


1.6x10* 


am 


ps/m2,  Oc=    12-16  mhos/m,  8-    .1m   (half  the  channel   width 


for   fully   developed  flow),  T      =2500<>K,  T     =18000,  The 

core  wall 

voltage  drop  parameter  is  0.282.  This  equates  to   a   voltage 

loss  per  boundary  layer  of  119  volts  forC7  =  16,  and  159  volts 
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Figure  A. 4  Non-dimensional  voltage  drop  4>   versus  wall  temperature  and 
free  stream  temperature  for  Cesium  seed,  constant  c  . 
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Figure  A.  5    Non-dimensional   voltage   drop    (J)     versus  wall    temperature   and 
free   stream  temperature   for  Cesium  seed,    varying  C    . 
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Figure   A. 6     Observed   transverse  voltage  distribution. 
(After   Sonju  and  Teno3°) 
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for(7t=12.   Extending  the  voltage  profile  to  the  channel  wall 

in  Fig.  A6  shows  a  V    of  about  150  volts. 

bl 


Another  example  comes  from  Argyropoulos  et  al[ 33  ]  who 
use  a  sophistocated  computer  technique  to  compare  theory 
with  experiments.  They  arrived  at  a  boundary  layer 
potential  drop  at  the  anode  of  101  volts  for  an  experiment 
on  the   AVCO-APL   channel.    Using   the   same   data,   namely 

T     =  2600OK,     T     =  2000OK,   Q   =  1.02,    (5=60  mm,    a 
core  wall  ^ 

chemistry  of  toluene,  oxygen  and  cesium,  and  the  1/7th-power 

law,   Fig.   A. 5   gives   a   non-dimensional  potential  drop  of 
0.24.   This  corresponds  to  a  predicted  drop  of  110.5   volts. 

The   next   example   shows   the  importance  of  selecting  a 

proper  model  for  the  boundary  layer.   Fig.  A7  is   from   data 

of   Kessler  and  Eustis[6].  The  parameters  for  the  experiment 

were:   I  =  2.75  amps,  electrode  surface   area  =   3.68   cm2, 

T     =  16850K,  T     =2700OK,  6c  =  10  mhos/m,  P=1.5,   The  gas 
elec  core 

was   Ethanol   in   oxygen   with   155   KOH   added   to   increase 

conductivity.   Nitrogen  was  introduced  for  cooling  such  that 

the  N  /O   ratio  was  0.5.   Using  the  1/7th-power  law  cf>    was 
2   2  bl 

found   to   be   0.290.    Kessler  predicted  the  boundary  layer 

thickness  to  be  1.2  cm.   The  calculated  voltage  drop  is  then 

8.5   volts  falling  well  below  the  46  volts  shown  in  Fig.  A7. 

The  difference   is  too  great   to   be   attributed   to   sheath 

phenomena  alcne. 

The  geometry  of  this  last  case  is  the  cause  of  the 
discrepency.  The  channel  depth  is  3.1  cm  while  the  width 
between  opposing  electrodes  is  10  cm.  This  means  that  the 
flow  does  not  behave  as  between  two  flat  plates  since  the 
sidewalls  have  a  strong  effect  on  the  electrode  wall 
boundary   layer.    The   model   used   in   this   comparison  is 
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(After  Kessler  and  Eustis   ) 
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obviously  inappropriate.  Schlichting  (Ref.  38,  p.  575) 
contains  information  on  the  problem  of  wall  interference  in 
non-circular  cross-sections  and  points  out  the  increase  in 
friction  and  boundary  layer  size.  A  gross,  but  closer 
approximation  for  this  procedure  is  to  assume  fully 
developed  flow  (as  was  suggested  by  Reseck[U3]  for  this  same 
channel).  This  results  in  a  voltage  drop  prediction  of  35.5 
volts.  Kessler  attributes  10  volts  to  the  sheath  drop  for 
the  cathode.  When  combined  with  the  35.5  volts  predicted 
above,  the  entire  drop  is  accounted  for  more  fully. 

8«   Conclusions 

Equation  A. 20  represents  a  useful  means  of  determining 
thermal  boundary  layer  voltage  losses  as  long  as  the 
following  conditions  are  met: 

1)  A  function  of  enthalpy  vs.  temperature   is   available, 

or,   for   the   case   of   Eg.    A. 16,   the  gas  has  a  fairly 

constant  C   in  the  range  of  temperatures  of  the  channel. 
F 

2)  the  core  conductivity  is  known  or  can  be  estimated  with 
some  degree  of  accuracy, 

3)  Hall   parameter,  core  temperature  and  wall  temperature 
are  known,  and 

4)  a  suitable  boundary  layer  description  is  available. 

The  simplicity  of  this  approach  makes  it  desirable  for 
integration  with  a  computer  program  which  predicts  the 
performance  of  an  HHD  generator.  Such  programs  usually 
already  contain  free  stream  and  wall  temperatures,  core 
conductivity,  and  Hall  parameters.  Since  these  eguations 
are  one-dimensional  they  allow  calculations  at  stations 
along  the  channel,  eliminating  the  need  to  determine 
properties  at  points  in  a  two-dimensional  field.  A  sample 
program  applying  this  method  at  one  station  is  given  in  the 
section  "Computer  Programs." 


109 


APPENDIX    B 

£EOOF    OF    THE    NON-EXISTENCE    OF    A    ONE-DT MENSIONA L    SHEATH 

SOLUTION 

A  first  attempt  at  solving  a  complicated  set  of 
equations  is  usually  to  reduce  the  problem  to  one  dimension. 
For  the  equations  which  describe  the  sheath  this  procedure 
seems  reasonable  due  to  the  small  thickness  of  the  sheath 
relative  to  the  electrode  dimensions.  Such  an  attempt  for  a 
constant  temperature  plasma  with  no  ionization  or 
recombination  will  fail  because  of  the  nonexistence  of  such 
a  solution.  Results  published  in  Rei.  21  are  summarized 
here  along  with  significant  conclusions  which  were 
instrumental  in  shaping  the  course  of  the  two-dimensional 
work. 


In   order  to  simplify  the  analysis  a  change  of  variables 
is  made. 


r  =  n.  -  n 
1    e 


£  =  n .  +  n 
s    1    e 


(B.1) 


The  following  manipulations  are  made  on  Eqs.  26   und   27 
the  text: 

1)  Magnetic  field  effects  are  dropped, 

2)  The  equations  are  considered  in  only  one-dimension, 

3)  The  equations  are  added  to  each  other  and  subtracted, 

4)  The  change  of  variables  in  Eg.  B. 1  are  applied. 
The  results,  including  Eg.  8  are 


in 


£-4  =  -  r 

2 

dy"  (   dy)  +  dy2 


=   0 


(B.2) 


(B.3) 
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S_  u   Si)  +  Sfr  =  o 

dy  (^  dy)    dy2      ° 


(B.4) 


If  Egs.  B.3  and  B.4  above  are  integrated  and  E  =  -  d4>/dy 
is  introduced  then 


dE 

dy      ' 


(B.5) 


rE  +  Si   =   C. 
dy       1 


(B.6) 


-  £E  + 


dr 

dy 


=   C. 


(B.7) 


It   is   not   difficult  to  show  that,  in  the  absence  of  a  net 
ion  current 

Cl  =  ~C2  =  V(eV  (B"6> 

where  J   is  the  conventional  current   flux   and   D    is   the 
Y  e 

electron  diffusion  coefficient.   Obviously,  since  J   exists, 

y 

C   and  C   are  non-zero  constants. 
1       2 


Because    the    ambipolar   region   is   quasi-neutral,  l 
approaches  zero  away  from  the  anode,   reaching   a   value   of 
zero  at  the  flasma  proper.   Thus,  as  y  increases: 

P  approaches  0 

F   approaches  F 

E  approaches  E^ 


Eventually  a  distance  y  =  y   would  be  reached   so 

o 


every  y  >  y 


| TE |  <  Cx/2 


that   for 


(B.9) 


Hence 
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§£  >  V2 


for  all   y  >  y 


(B. 10) 


S(y)  >  C(yQ)  +  (c^/2)  (y-y  ) 


(B-11) 


which   does    not    approach 
Furthermore,  since 


constant   as    required. 


dy     ul    ^ 


(B. 12) 


g  >  -Cx  +  [?(yo)  +  Cl/2  x  (y-yo)]  E 


(B.13) 


for  y  >  y  /  which  eventually  grows  at  least  as  fast  as 
c 

C 

(B. 14) 


»-£  (y-y0) 


and  does  not  approach  zero  as  required. 


These   cbvious   inconsistencies   arise   because   of   the 

constant  C  .   The  necessary  requirement  for   a   solution   is 
1 

that  the  current  density  decrease  away  from  the  anode  giving 

rise  to  current  constrictions  at  the  surface.    Under   these 

conditions,  C   would  decrease  at  an  appropriate  rate  so  that 
1 

P  approaches  zero  and  c  approaches  F    .    Mechanisms   which 

allow  a  proper  variation  of  C   are  diffusion  in  two  or  three 

dimensions  and  ionization/recombination.  The  former  negates 
the  applicability  of  a  one-dimensional  Cartesian  solution 
while  the  latter  eliminates  the  frozen  flow  assumption. 
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APPENDIX  C 

DIMENSIONAL  AND  FRACTIONAL  ANALYSES  OF  THE  SHEATH  EQUATIONS 

Equations  1-6  in  the  text  are  cumbersome  in  that  they 
contain  many  physical  parameters  which  in  their  present  for  a 
must  be  carried  along  as  the  equations  are  operated  on.  It 
is  useful  to  conduct  a  dimensional  analysis  or  these 
equations  in  an  attempt  to  collect  the  coefficients  and 
variables  in  such  a  way  as  to  simplify  the  appearance  of  the 
equations  and  better  understand  their  nature.  Also,  before 
introducing  the  problem  into  the  computer  it  is  imperative 
to  know  the  sizes  of  the  quantities  that  are  being  dealt 
with.  For  this  latter  reason,  a  fractional  analysis  showing 
the  relative  significance  of  certain  physical  effects  is 
also  given  in  this  appendix. 

1 •   Di^ensio nal  Analysis 

The  method  outlined  by  Langhaar[ 44  ]  is  utilized  here   to 

determine   the   non-dimensional   parameters   relevant  to  the 

sheath  equations.   The  significant  parameters   which   appear 

in   Eqs.   1   and  3  are  J,  <j>,  n  ,    y,  e,  £  ,  D  ,  U  ,  Q ,    and  B, 

s         o   s   s  ^ 

where  y  represents  any  coordinate  length.   An  expansion   and 

manipulation   of   Eg.  1  with  Eqs.  3  and  5  will  show  that  the 


Since   k   and  T 


o 


terms  D   and  U     can  be  replaced  by   kT  . 
s      s  o 

always   appear   together,   they  can  be  expressed  as  a  single 

term.    Additionally,   the   Hall   parameter  ft      is    already 

non-dimensional   and   will  not  be  included  in  this  analysis. 

This  leaves  <D,  n  ,  y,  e,  £  ,  and  kT  .    In   the   MKS   system 

s         o        o 

these   parameters  are  combinations  of  the  units  of  -mass  (LI)  , 

length (L),  charge  (g)  and  time  (t) .   The  matrix  showing   the 
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coefficients  of  these  units  for  each  parameter  is 

<t>     n     L     e    £     kl0 

o 

0     0 

P 

0     0 


M 
I 


1 

2 
-  1 


11 


-2 


-3 

0 


0  j  -3 


2 

2 


1  I 
2 

0 
2 


L 


At   first   glance   it   might   be   expected   that   there   are 

6-4=2   dimensionless   parameters   (number   of  parameters 

minus  the  number   of   units   necessary   to   describe   them) . 

However,   closer   examination   shows   that   the  t  row  is  not 

linearly  independent  of  the  M  row.   In  fact  the  rank  of  this 

matrix    is   three   not   four.    Thus   there   are   6-3=3 

independent  dimensionless  parameters  (number   of   parameters 

minus   the   rank  of  the  matrix)  which  describe  the  system  of 

equations.   This  result  would  also  have  arisen  if  k   and   T 

o 

had  been  considered  separately. 

Disregarding  the  t  row,  the  homogeneous  linear  algebraic 
equations  whose  coefficients  are  the  numbers  in  the  rows  of 
the  dimensicr.al  matrix  are 


al  "   a5  +  a6   =   ° 
2a1  -  3a2  +  a3  -  3a5  +  2ag   =   0 

-a1  +   a4  +  2a5   =   0 


(C.1) 


Solving  for  a  ,    a  .  and  a   in  terms  of  the  other  three  gives 
4    5       6 
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a.  =   a-,  +  6a~  -  2a., 


a5  =  -3a2  +  a3 


(C.2) 


a  —  —   3-i  ™*  J3n  t  a~ 
6      ±23 


Now  setting  the  values  a  =1,  and   a  =a  =0,   the   result   i: 

1  2   3 

a  =1,  a  =0,  and  a  =-1. 

4      5         6 


Continuing  by  assigning  a  value  of  one   successively   to 

a    and  a   with  the  other  two  of  that  set  equal  to  zero,  the 
2        3 


solution  matrix  is  generated. 


+ 


4>    n     L    e 

_U i 


kOL 


1 


0  !      0 


j 


I   i 

i 


... 


-1 


3 


ILL 

Li  ° 


i 


0 


1  I  -2 

_L 


The    resulting    non-dimensional    parameters   are    then 


e< 
kT 


o 


,   e2  ,3 

n  =  (kT-T-)  n 
o  o 


e  kT 
o   o 


y 


(C.3) 


*       A 


These   terms   are  all  independent  since  $,    n  and  y  appear  in 
only  one  of  each.   In  numerical  form  they  look  like 


J  =  $/<{,   =  I.16xl04  T   1  4) 
n  =  n/n   =  9.26xl0~12  TQ"3  n 


(C4) 


y  =  y/L  =  4. 76x10 J  TQ  y 
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2.   Estimate  of  the  Sheath  Size 

The   sheath   size  can  be  estimated  by  a  technique  called 
fractional  analysis  or  approximation  theory[45].   Sq.  3  is 


V  *  =  "(n.  _ne)e/eo 


(C.5) 


Now  "characteristic"  parameters  <J>  ,  \    ,  and  n   are   factored 

os       o 

out  of  the  equation  for  each  variable. 


(C-6) 


where   <J>    and  n   are  the  maximum  value  of  the  potential  and 


o 


charge   concentration   respectively,   and  \       is   a   sheath 

s 

charactersitic        length.  Since      in      the      sheath      all      the 

variables   are       of      order      unity,       for      the      equation      to      be 

n0  eA 

consistent   tne   term   -r must   also   bo  of  order  unitv. 

<£o 

Setting  this  term  equal  to   one   and   solving   for  A  ,   the 


result  is 


X. 


t   Yo  o 
noe 


(C.7) 


This  gives  an  order  of  magnitude  for  the  characteristic 
sheath  length.  For  example,  for  a  sheath  drop  of  10  volts 
in  a  plasma  ionized  to  1Ci8  electrons/m3  the  characteristic 
sheath  length  is  2.4x10-sm. 

3.   Characteristic  Current  from  Conduction 


From  Eg.  10  the  conduction  current  is 


J   =  -  u  en  V(f> 
e     e   e 


(C.8) 
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Working  in  one-dimension  and   identifying   a   characteristic 


current  J  ,  a  fractional  analysis  gives  the  equation 
o 


J  X 
o  s 


<b    u  en 
Ycv  e   o 


*  1 


(C.9) 


Notice   that   in   order   to  be  consistent  the  characteristic 

length  in  the  potential  derivative  must   again   be  A  ,   the 

s 

sheath   length.    Solving   for   J   when  n   =  10iam-3/  (J     =  1 

o       o  e 

in2  (volt-sec)  -1,    and    <J>   =    10,    gives    J    =    6.8x10*   amps/m2.         The 

o  o 

value   for  fj.     comes  from  Eg.  16  for  electrons  colliding  with 
e 

CO   molecules  at  1000OK. 
2 

**•   Diffusion  L§IL3£ii 

The  potential  is  a  monotonic  positive  function  which  is 
suitable  for  approximation  theory,  however  the  charge 
density  profiles  have  many  inflection  points  in  the  sheath 
and  are  not  well  suited.  Nevertheless,  it  is  possible  to 
consider  profiles  beyond  the  sheath  where  characteristic 
lengths  can  fce  defined. 

If  diffusion  is  to  play  a  significant  role  in  current 
production  then 


J  - u  kT  |£ 

o   He   e  3y 


Again,     using    fractional 
approximation  is  made: 


o  a 


(C.10) 
analysis,    the    following 

(C.  11) 


n    uJ<T 
cr  e      e 


«  1 


Solving    for    I      with    the    same    parameters    as    above   gives   L        = 
d  a 

2.03x10-7.         Since      L   <A       it    can    be    concluded   that    any    role 

d      s 

played  by  the   mechanism   of   diffusion   should   be   visible 
within  the  sheath. 
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5-   Ionizatio^Recombinjition  Length 


From     Hinnov    and    Kirschberg[  46  ]    the    three-bo<iy 
recombination  rate  of  electrons  and  ions  is  given  by 

dn_ 


dt 


a  n 


(C.12) 


where 


c  c    -,n-33  ,kT.-4.5 
a  =  5. 6x10     ( — )      n 

e        e 


otic 


(C.13) 


Using  the  technique  outlined  previously, 

no  dne    c  a    ln-33  ,kT\-4.5    3-3 
r~  7T-  =  5.6x10     ( — )      n    n 
x   dt  e        o    e 


(C.14) 


where  T     is  the  characteristic  recombination  time.   Let   the 
r 

term   5. 6x1 0~3<)  (kT/e)  -♦*  snz   be  of   the   order  of  unity  aa.d 


solve  for  T  . 
r 


_-3->    VT-45    ?  -1 
Tr  -  (5.6x10  JJ)  {£*.)  (nQZ) 


(C.15) 


At   T=1000OK   and   n   =  1Ql8m-3,   the   characteristic    time 

° 

becomes  2.9x10-3sec. 


In  order  to  change  the  characteristic  time  to  a 
characteristic  length  the  mean  thermal  velocity  is  used 
since  this  is  the  transport  mechanism  for  ionization  and 
recombination.   Thus 


L   =  t  c  =  T 
r    r     r 


(C. 16) 


Tim 


Using   the   same   conditions  as  above  L   =  5.7x102  m.   Since 

r 

L    >>\      it   can   be    concluded    that    ionization    and   recombination 
r        s 
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do  not  have  time  to  occur  within  the  sheath  for  the 
conditions  used  in  this  analysis.  Since  the  results  of  this 
work  show  that  the  ambipolar  region  is  of  the  same  order  as 
the  sheath  length,  ionization/recombination  is  not  expected 
to  play  a  part  in  that  region  either. 

6.   Convection  Time 


At  atmospheric  pressure  and  temperatures  around  10000K 
sonic  velocity  is  of  the  order  of  103  ia/sec.  This  is  the 
order  of  velocity  of  the  flows  of  MHD  generators.  A  rough 
characteristic  time  relevant  to  the  sheath  could  ue 
established  by  comparing  the  flow  velocity  with  the  size  of 
the  sheath. 


t   =  X   /U 
c    s   c 


(C17) 


Using    4>      =    10   volts    and    n      =    KM8    m"3,    J     =    2.4x1G-8sec. 
o  o  c 


A    characteristic    sheath    time    can    be   defined    from 


s  s' 


(C. 18) 


for   comparison  with  T  .   Using  representative  values  for  \ 

c  s 

and  c,     T    =    1.2x10-*°  sec.   Thus,  any  effects  occuring   from 
s 

convective   phenomena   such   as   mixing   and  boundary  layer 
effects  do  not  have  sufficient  time  to  affect  the  sheath. 
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COMPUTER  PROGRAMS 


The  principal  program  used  in  this  work  is  the  one  that 
calculates  the  potential  and  charge  density  distributions 
for  the  sheath  and  ambipolar  regions.  This  program  is 
explained  in  Part  1  of  this  section  and  reproduced  there. 

Auxiliary  programs  are  those  that  calculate  the  current 
distribution  and  produce  the  outputs  used  in  this  work. 
Part  2  shows  the  auxiliary  programs. 

Part  3  is  a  reproduction  of  the  program  usei  in  the 
calculation  of  the  boundary  layer  in  Appendix  A.  A  sample 
enthalpy  table  is  provided  in  data  format  for  the 
Toluene/oxygen/cesium  mixture. 

1 •   Sheath  Program 

The  sheath  program  begins  with  an  assumed  solution  and 
operates  on  the  field  of  parameters  to  replace  the  assumed 
solution  with  a  more  accurate  one.  Theoretically,  the  more 
iterations  that  are  used,  the  more  accurate  the  result  is. 
An  example  of  an  initial  assumed  solution  is  P{I,J)=0, 
G(I,J)=BZ,  SI(I,J)=BZ,  and  TH(I,J)=1.0  throughout  the  field, 
(see  the  program  comment  cards  for  definitions)  These  three 
arrays  are  operated  on  and  a  more  accurate  solution  is 
stored  for  later  analysis  or  further  iterations. 

The  following  pages  are  a  reprint  of  this  program,  and 
an  explanation  of  the  parameters  and  subroutines  is 
included . 
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c 
c 
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SHEATH  COMPUTATION  PROGRAM 
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CALL    LCAD 
2    CONTINUE 
OUTPUT    AND    STORE    RESULTS 
WRITE(2)     P 
WRITE(2)    G 
WRITE(2)     SI 
WRITE(6,100  ) 
100    F0PMAT(3X, • DATA    HAS    BEEN    READ    INTO    FILE') 

200  FCRMAT(6D10.0) 
300    F0PMAT(4I4) 

201  FGFMATOX, 6012.3,414) 
STCP 

END 


SUER 
IMFL 
REAL 
COMM 
COMM 
CONM 
CO  MM 

COMM 
COMM 
DIME 
EQUI 
HS  =  1 
HR  =  . 
H4  =  H 
BETR 
NCCW 
NFLA 
NCCL 
NRCrt 
MUD  = 
MLD= 

BEGIN  I 
DO  9 
ND  =  2 

ROTATIO 
IN=N 
IF(  I 
IF(G 
CALL 

SWEEP  I 

30  DC  1 
I  =  IX 
I  F  {  I 
DC  9 
DO  9 
A(L, 

90  CONT 
I  F  ( I 
DO  2 

INTROCU 
TT  =  T 
THS  = 

LOAD  CC 
DO  2 
M  =  3* 
GO  T 
IF(J 
IF(  J 
IF(  I 
I  F  (  I 
IF(J 
DTY  = 
TS  =  . 
GO  T 

51  I F(  I 


LT 

CI 

4 

N/ 

N/ 

N/ 

N/ 

1) 

N/ 

N/ 

SI 

£L 

CO 

D 

/4 

BE 

NC 

=  1 

3  = 

3  = 


INE  LOAD 

T  REAL*t ( A-H,0-Z) 

EPS 

BH/P<  51,51)  ,G(51 

BL2/BET,GAM,H ,PHI 
3L3/BMAX»W, JMIN, JMAX 
3L4/FP(51),PPP(51),GXP(51) 


,51), SI (5 

, 6Z,NCGLS 


1,51) 

,NRGWS,NCFN,MAX 

,GXPP(51),SP(51) ,SPP( 


BL5/NLFT,NRIT 
BL5/TH(51 ,51) 
ON  T(153) , 
ENCE  (A(l, 

00/(H*H) 
OC/H 
.D  00 

T*HR/2.D  00 
OLS-1 

10 
NR0WS*3 


DTX(51) 

A(10,153)  ,B( 1505) 
1),B(1)  ) 


ERATIONS 
JCNT=1,.MAX 
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N  PARAMETERS 
FLAG-iNFLAG/4)*4 

N.EQ.0.CR.IM.EG.3)  N0=1 
AM.LT.l .0  00)  GO  TO  30 

TMP(ND) 
N  I  DIRECTION 

IX=1,NR0WS 

N.EQ.0.GR.IN.EQ.2)  I=NR0WS-IX+1 

0    M=1,NRCWS3 

0    L=1,NCCLS3 

M)=0.0 

INUE 

•  NE.l .AND.I.NE.MROWS)     CALL    DERV(I,2) 

J=1,NC0LS 
CE    TEMPERATURE     PARAMETERS 
H(  I,  J) 
TT:rTT*H  S 
MFUTATICNAL    MATRIX 

K  =  l,3 

(  j-i)+k 

0    (50,5  1),ND 

.EQ.l.AND.I.EQ.NCEN)    GOTO     (76, 21,21), K 

.EC.NCOLS)    GO    TO     (21, 19, 19), K 

.FQ.l)     GO    TO    91 

.EQ.NROWS)     GO    TO     (96,  97, 98), K 

.EG.l )     GO    TO    (40,19,19) , K 

TH( I , J*l) -TH( I , J-l  ) 

5D    00*TT*DTY*H4-BET*TT*HR*(PP( J)+.5D    00*DTX(J)) 

0     (5,6,7) ,K 

.EQ.l.AND.J. EQ.NCEN)    GO    TO    (76,  21,  21), K 
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INT 
5 


I  F  ( I  . 

IF(J. 
IF(  J. 
IF(I  . 
DTY=T 
TS  =  .5 
GO  TO 
ERNAL 
A(7,M 
A  (  8  ,  M 
A  (  3 ,  M 
A(9,  M 
A  {  6 ,  M 
T(M)  = 
GO  TO 
A  (  3 ,  M 
A(6,M 
A(9, 
A(  2, 


A(8 


T(l*)  = 


60U 
76 

91 


92 
93 
94 
96 
97 
98 
86 
87 
88 
89 
21 
19 
40 


T 
GO  TO 
A  (  3 ,  M 
A  (  6  t  M 
A(9t  M 
A  (  1 1  M 
A  (  7  ,  M 
T(M)  = 
GO  TO 
NDARY 
T  (  M )  = 
GO  TO 
I  1=1 
J1  =  J 
IF(ND 
IF(ND 
GO  TO 
T(N>  = 
GO  TO 
T(M)  = 
GC  TO 
T  (  M )  = 
GO  TO 
T  (  M )  = 
GO  TO 
T(P)  = 
GO  TO 
T  {  M )  = 
GC  TO 
A  (  3  »  M 
GO  TO 
T(M)  = 
GO  TO 
T(M)  = 
GO  TO 
T  (  M )  = 
GO  TO 
T(M)  = 
GC  TO 
T  (  M  )  = 
GO  TO 
A(6,M 
A(9,  M 
A<  7,M 
At  10t 
T(M)  = 
GO  TO 


EQ.:\!RO 
Et.l  ) 
EQ.NCG 
EQ.i) 
h(I , J+ 
C  OG*T 

(5,6, 

FIELD 
)=-l.D 
)=1.D 
)=HS 
)=HS 
J =-4.0 
-PPP( J 

2 
)=THS- 
)=-4.D 
)=THS+ 
)=TT*( 
ET#HR* 
)=-A(2 
TT*(-T 
5D  00* 
5D  00* 
5D  0  0* 
T*GXP( 

2 
)=HS 
)=-4.D 
}=HS 
)=-( SI 
)=-<Ul 
SI(  I  ,J 

2 

CONDI 
PHI 

85 


WS) 

GO 

IS) 

GO 

1)- 

T*D 

7), 
CO 
00 

00 


GO  TO  1 
TO  91 

GO    TO    86 
TO    (4i,19,19),K 
TH( I, J-l ) 
TY*H4 
K 
EFFICIENTS 


00*HS 
) 

TS 

C0*THS 
TS 
G( 
<T 
,M 
T* 
CT 
G( 
BE 
J) 


I  ,  J* 
T*GX 
) 

GXPP 
X(  J  ) 
I,  J) 
T*OT 
) 


1)-G( I , J-l)  )*H4-.5D  00*G( I , J)*DTY*H4- 
P(J)-.5D  00*G( If J)*DTX(J)J 

(J) +GXP( J)*PP< J)- 

*GXP< J)+G< If J)*(G(  I  ,J)-SI< I ,J)  ))- 

*DTX( J)*PP( J)- 

Y*HR*(G( I  ,  J)*PP< J)  + 


00#HS 

(  I,J  +  1)-SI( I ,  J-l)  )*H4 

,  M) 

)*(SIU,J)-G(I,J))-PP(J)*SPU)-SPP(J) 

TIONS 


.EQ.I)  I1=NR0WS 
.EC. 2)  J1=NC0LS 

(9  2,9  3,94) ,K 
P( II, Jl) 

85 
G(Ilf Jl) 

p5 
SI( II, Jl) 

85 
P(2V J)-P(lt J)+P(NCOM, J) 

85 
G(2,  J)-G(  1,  J)+G(NCOM,  J) 

85 
SI(2,J  )-SI( 1 , J  J  +SI (NCOM, J) 

85 
)=-l.D  00 

(87,  88, 89), K 
P(I  ,2)-P(  If  1) 

85 
G(  I,2)-G( 1,1) 

85 
SI(I,2)-SI(  I,  1) 

85 
CD 

85 
EZ 

85 
J-3.D  00*BZ 
)=-4.D  00*8Z 
)=-3.0  00 
M=4«D  00 
-BZ*P(  I,3)+G(  I,3)+BET*BZ*(  P( 1  + 1  ,  1  )-P( I  -  1 ,  1  )  ) 

2 
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200 


202 


41  T(M) 

85  A(6, 

2  CQNT 

MODIFY 
KMAX 
MB  =  0 
MK=6 
DO  2 
MK  =  M 
IF(K 
DC  2 
IF(K 
M=10 
MA  =  M 
6  (MA 
COM 
KMM= 
MK  =  0 
DO  2 
MB=M 
MK=M 
MP  =  1 
DO  2 
IF(K 
M=10 
MA  =  M 
B(PA 
COM 

SOLVE  M 
CALL 
IP*I 
IF(  I 
IF(  I 

SUCCFSS 
DO  1 
M  =  2* 
HINT 
HIM 
HINT 
P(  I, 
G(  I  , 
SI  (I 
CO  NT 

PRINTER 
WRIT 
IF(J 
WRIT 
CALL 
WRIT 
CALL 
WRIT 
CALL 
IF(  I 
CALL 
CALL 
CALL 
CALL 
NFLA 
RETU 


=  P(  I, J) 

M)=1.D  00 

INUF 

MATRIX  FOR 
=NPCWS3-4 


USE  IN  SUBROUTINE  DGELB 


00  KA=1 ,KMAX 

K-l 

A.LT.6)  MB=MS+MK 

00  KB=1,10 

A.LT.6. AND. KB. LE.MK)  GO  TO  200 

*<KA-1) +KB 

-MB 

)=A(KBtKA) 

INUE 

KMAX+1 


13 
1 


43 


99 


100 
101 
102 
103 
104 


FORM 
FORM 
FORM 
FORM 
FORM 
END 


02 
K  + 
K  + 
0- 
02 
B. 
*( 
-M 
)  = 
IN 
AT 
C 
+  1 
N  . 
P. 
IV 
3 
(J 

1  = 

2  = 

3  = 
J) 
J) 

,J 

IK 

C 

E( 

CN 

E( 

W 

E( 

W 

E( 

W 

N. 

P 

P 

R 

R 

G  = 

RN 

AT 

AT 

AT 

AT 

AT 


KA=KFM,NR0WS3 
MB 
1 
MK 

KB=1,10 
GT.MPJ  GO  TO  202 
KA-1) fKB 
B 

A(K3,KA) 
UE 

(JTV 

GELB( T,B,NR0WS3, 1 , MUD,MLO  ,  EPS ,  IER) 


EQ.O. 

GT.l  . 

E  OVE 

J=1,N 

-1J 

T(M+1 

T(M*2 

T(M  +  3 

=  P(  I, 

=  G(  1  , 

)=SI( 

UE 

UTPUT 

6,103 

T.LT. 

6,100 

PITER 

6,101 

RITER 

6,102 

RITER 

EQ*  1. 

CTAT( 

CTAT( 

CTAT< 

CTAT( 

NFLAG 


OR.IN.EQ. 
AND.IP.LT 
R-RELAXAT 
COLS 

}-P(  I  ,  J) 
)-G(I , J) 
)-SI( I, J) 
J) +W*HINT 
J) +W*HINT 
I, J)+W*HI 


2)  IP=I-1 

.NROWS)  CALL  DERV(IP,1) 

IGN  FEATURE 


1 
2 

NT  3 


)  JCN 

J  M  IN. 
)  JCN 
(  P,NR 
)  JCN 

(CNR 


)  JCN 
(SI,  N 
OR.  IN 
P,NRO 
G.NKO 
SI  ,NR 
TH,NR 
+  1 


T 
OR.  J 

T 

OW  S , 

T 

OWS, 

T 

*OWS 

,EQ. 

WS,N 

WS,N 

OWS, 

OWS, 


CNT.GT.JMAX)  GO  TO  43 

NCCLS) 

NCOLS ) 

,NCOLS) 

3)  GO  TO  99 

COLS) 

COLS) 

NCOLS) 

NCOLS) 


(//,3X, 'POTENTIAL  PLOT*, 14,/) 
' ON  PLOT' ,14,/) 
0T«  ,  14,/) 


(  /  /  ,  3  X  , 
(//,3X, 
(3X, 14) 
(//,3X, 


•  ELECTR! 
•ION    PL  i 


•TEMPERATURE     PROF  I LE«  , 14  ,/ ) 
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SUBROUTINE  DERV( I, NF ) 

IMPLICIT  *EAL*8( A-H,0-Z) 

COMMCN/BL1/ P< 51, 51 ) ,G( 51,51),  SI (51,51) 

CGMMON/BL2/ BET , GAM, H , PHI , 

C0MM0N/BL4/PP(51 I ,PPP(51I 


51) ,DTX(51 ) 


51) 
COMMGN/BL6/TH(51, 
LM=NROhS-l 
11=1+] 
12=1-1 

IF(I.EQ.l)  I2=LM 
IF( I .EQ.NROWS)  11=2 
JF(NF.EG.2  )  GO  TO  4 
FIRST  DERIVATIVES 
H2=.5D  OO/H 
DC  1  J=l,NCOLS 

DTX<  J)=(TH<  II,  J)-TH(  12,  J)  )*H2 
PP(J)=(P(I1,J)-P(I2,J) )*H2 
GXP(J)=(G(I1,J)-G(I2,J) )*H2 

1  SP(J)  =  (SI(  II,  J)-SI(  12, J)  )*H2 
RETURN 

SECOND    DERIVATIVES 
4    HS=1 .CD    00/(H*H) 
DC    2     J=1,NCCLS 
P  P  P  (  J  )  =  (  P  (  1  1  ,  J  )  +  P  (  1 2  ,  J  )  K-  HS 
GXPP( J)=(G( II, J5+G( 12, J) )*HS 

2  SPP( JJ=(SI( I1,J)+SI( 12, J) )*HS 
RETURN 

END 


BZ,NCOLS,NRCWS,NCEN,MAX 
,GXP( 51) ,GXPP(51) , SP( 51) ,SPP( 


SUBROUTINE  WRITER! Z , NCOLS,NROWS ) 

IMPLICIT    REAL*8< A-H,0-Z) 

DIMENSION    Z(NROWS,NCOLS) 

J  1  =  1 

J2  =  l  1 
10    IF(J2.GE.NCCLS)     J2=NC0LS 

DO    15    I=1,NR0WS 
15    WRITE(6,220)     I ,  ( Z ( I , J )  ,  J  =  J 1 , J2 ) 

WR1TE(6,240) 

IF(  J2.ECNC0LS)     RETURN 

J1=J1+11 

J2=J2+11 

GO  TO  10 
22  0  FORM  AT (•  »,I3,11EU.4> 
240  FORMAT!/) 

END 
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SUBROUTINE  ROT  AT ( A , NROWS » NCOLS ) 

IMPLICIT  REALMS  (A-H,0-Z) 

DIMENSION  A (NROWS. NCOLS) ,B( 51,51) 

DO  1  1=1, NROWS 

DO  1  J=l, NCOLS 

B( I  * J)=A(  I, J) 

DO  2  1=1, NROWS 

DO  2  J=l, NCOLS 

A< I, J)=B( J,  I) 

RETURN 

END 


SU 
IM 

CO 

cc 

CG 
NR 
NC 

Gl 

G2 
DO 
DO 
EL 
DP 
DP 
DN 
ON 
AR 
AR 
TH 
CO 
IF 
TF 

DO 
TH 
Th 
DO 
Tri 
TH 
RE 
TF 

DO 
TH 
TH 
DO 
TH 
TH 
RE 
EN 


eROUT 
PLICI 
MMCN/ 
KiiON/ 
PMCN/ 
CM  =  NP 
CM=NC 
=GAM* 
=GAN/ 
2  1  = 

2  J  = 
=G(  It 
X=P(I 
Y=P(I 
X=G(  I 
Y  =  G(I 
G=G1* 
G2  =  G2 
(I  ,  J) 
NTINU 
(  N  C  .  E 
=  (  1  .D 

•  0 

3  J  = 
(NROW 
(1,  J) 

4  1  = 
(I  ,NC 
(1,1) 
TURN 

=  (  i.D 
.D 

6  1  = 
(1,1) 
{ I  ,  NC 

7  J  = 
(NRCW 
(1  ,  J) 
TURN 
D 


INE 
T  RE 
BL1/ 
BL2/ 

BL6/ 
OWS- 
GLS- 
.250 
(H*H 
2,NR 
2,NC 
J) 

+  1,J 
,  J  +  l 
+  1  ,  J 
t  J*l 
(  DNX 
*(  DP 
=  (-A 

Q.l  ) 

00  + 
00 

1,NC 
S,  J) 
=  1  .D 
2,NR 
OLS) 
=  TH( 


(NO) 

8( A-H,0-Z) 

1,51) ,G( 51 ,51 ) ,SI (51,51) 

,GAM,H, PH I, BZ,NCOLS,NROWS,NCEN, MAX 

51,51) ,0TX(51 > 


TMP 
AL* 
P(5 

BET 
TH( 
1 
1 

00/(H*H) 
) 

CM 
CM 

)-P( 1-1, J) 

)-P( I ,J-1 ) 

)-G(I-l,J) 

)-G( I , J-l ) 

*DPX+DNY*DPY)/EL-1.D    00 

X*DPX+DPY*DPY) 

RG+DSQRT( ARG-ARG+ARG2) )/2 


D   00 


GO    TO    5 
DSQP.TU  .D 


00 +4  .0    00*GAM*P(NRGM,NCEN)**2) )/2 


OLS 
=  TF 
00 
OM 

=  TH(I  ,1) 
I,2)-TH(  I,NC0LS)+TH( I,NCOM) 


C0+DSQRT(1 
00 

1 ,NRCWS 
=1.0  00 
0LS)=TF 
2,NCCM 
S, J)=TH(1  ,J 


D  00+4. D  00*GAM*P(NCEN,NC0M)**2))/2 


=  TH(  2,J)-TH( NROWS, J)+TH(NROM, J) 
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2.   Auxiliary  Programs 

Several  auxiliary  programs  were  used  to  process  the 
results  of  the  sheath  program.  The  first  calculates  the 
non-dimensicnal  current  from  the  potential,  charge  density 
and  temperature  profiles.  The  second  program  creates  a 
graphical  display  of  the  potential  and  charge  distributions 
along  a  cut  perpendicular  to  the  electrode  and  through  it. 
The  third  program  uses  the  same  cut  to  regenerate  the 
electron  temperature  in  a  one-dimensional  plot.  The  fourth 
program  is  a  modification  of  one  by  Biblarz[47]  which 
produces  a  two-dimensional  display  or  contour  map  of  the 
potential  and  charge  density  fields. 

These  programs  are  reproduced  in  the  following  pages  and 
include  explanations  of  their  parameters  and  use. 
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C**  CURRENT    PROGRAM  ** 

C 

c 

C   THIS  PPCGPAM  IS  CESIGNED  TO  CALCULATE  THE  CURRENT 

C   DISTRIBUTION  GIVEN  THE  POTENTIAL  AND  CHARGE  DENSITY 

C   FIELDS.   IT  OUTPUTS  USING  SUBROUTINE  WRITER  FROM  TFE 

C   SHEATH  PROGRAM.   AS  WITH  THE  SHEATH  PROGRAM  IT  USES  A 

C   DATA  CELL  FOR  VARIABLE  STORAGE.   IT  REQUIRES  THE 

C   FCLLCWING  PREDETERMINED  PARAMETERS: 

C      P      POTENTIAL  DISTRIBUTION 

C      G      ELECTRON  DENSITY  DISTRIBUTION 

C      SI     ION  DENSITY  DISTRIBUTION 

C      GAM    JOULE  HEATING  PARAMETER 

C      BET    HALL  PARAMETER 

C      ZL     STEP  SIZE  MULTIPLE 

C      NROWS  ARRAY  WIDTH 

C      NCOLS  ARRAY  LENGTH 

IMPLICIT  REAL*3( A-H,0-Z) 

DIMENSION  P(51,51),G<51t5l)fSI(51f51)fXJ<51,51)tYJ<5i, 
51)»SX(51,51),ZI(51),TH(51,51I 

READ(l)  P 

RE£D( 1)  G 

READ(l)  SI 

GAM=0.0 

ZL=5.C  02 

BET=1.D  CO 

NRGWS=51 

NCCLS=51 

NRCM  =  i\RCWS-l 

NCCM=,\C0LS-1 

NB=NCCLS~1 

NC=NE-1 

H=ZL/DFLCAT(NR0WS-1) 

HR=5.0D-01/H 

Gl=GAf-*.2  5D    00/(H*H) 

G2=GAf/(H*H) 

DO    10    1=1 i NROWS 

DC    10    J=l, NCOLS 
10    Th( I , J)=1.D    00 

IF(GAP.LE.l.D-04)    GO    TO    20 
C       CALCULATE    DERIVATIVES    AND    ELECTRON    TEMPERATURE 

DO    2    1=2, N ROM 

DC    2    J  =  2,NCCi-' 

EL=G(  I  ,J) 

DPX=F(I*1,  J)-P(  I-lt  J) 

DPY=P(I»  J  +  l  )-PU,J-l  ) 

DNX  =  G<H-1  ,J)-G(  I-lt  J  ) 

DNY=G( It J+l )-G(I t J-I ) 

ARG=G1*( CNX^DPX+DNY*DPY)/EL-1.D    00 

ARG2=G2*(  DFX*DPX+DPY*DPY) 

TH(I , J)=(-ARG+DSQRT( ARG*ARG+ARG2 ) )/2.D    00 

C       CALCULATE    CROSS-CHANNEL    CURRENT    DENSITY 

20    SI(IflTi-3?*GlItl)+4.*G(I|2)-G(I,3)-G(Ifl)*(- 
3.*P( I t 1)+4.*P<  I ,2)-P(  1.3)  ) 
SKI  , NCOLS )=G(  I,NC)-4.*G<  I,NB)+3.*G(  I,  NCOLS  )- 
GU,NC0LS)*(P(I»NC)-4.*P(I»NB)  + 
3.*P< I, NCOLS)) 

3  SI<I,J)=iH?I,J)*(G{I,J+l)-G(I,J-ll )-G(I,J)*(P(IiJ+l)- 

P(  I  f J-l  )  ) 

C       CALCULATE_AXIAL    CURRENT    DENSITY 

SX(l,jT=13"batJ)+4.*G(2,J)-G(3tJJ-G(l,J)*(- 

3.*P( 1» J)+4.*P< 2, J)-P(3,J) ) 
SX(NPCWS, J)=G(NC,J)-4.*G(NBt J ) +3. *G(NROWS » J )- 

G  (NROWS  T  J  )  *<  P(  NC  ,  J  )-4 .*P(  NB ,  J  )  ♦ 

3.*P(NRGWS» J) ) 

4  SX(I,J)=fH(I,J)*(G(H-l.J)-G<l*l,Jn-G{I,JI*(P(I*ltJ)- 

P(I-l.J) J 
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APPLY  H 
DC  5 
DC  5 
YJ(  I 
XJ(  I 

5  CONT 
INTEGRA 
FAC  = 
DO  7 
SUM1 
SUM2 
SUM= 
DO  8 
SUM1 
DO  9 

9  SUM? 
II  (I 

7  CCNT 
OUTPUT 
WRIT 
CALL 
WRIT 
CALL 
WRIT 
STCP 
FORM 
FORM 
FORM 
END 


ALL  PARAMETER 

I=l,NROWS 

J=1,NCCLS 
i J)  =  (SI (  I,  J)-BET*SX(  I  ,J) )*HR 
,  J)=( SX(  It J)+BET*SI( I t J) )*HR 
INUE 

TE  FOR  NET  CURRENT 
l.D  00/(3. D  00*DFLOAT(NB) ) 

I-1,NR0WS 
=C.O 
=  C.O 
XJ(I, 1)+XJ(I ,NCOLS) 


=SUM1*XJ(  I  ,  J  ) 

J=3,NC  ,2 
=SUM2+XJ( I TJ) 
)  =  (SUM+4.D  00* SUM  1  +  2. D 
INUE 


00*SUM2)*FAC 


101 
102 
104 


E(6t 101  ) 

WRITERU  J  ,  NRCWS  »NCOLS  ) 
E(6,102) 

WRITER (YJtNROWStNCOLS) 
E(6,104)  <  I,ZIU  )  ,I  =  1,NR0WS) 

AT(//,3X, 'X-CURRENT' ,/) 
AT(//f3X,«Y- CURRENT1 ,/ ) 
AT(//,51(3XtI4tiPD12.5f/) ) 
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c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


CNE-DIMENSIONAL  POTENTIAL  AND  CHARGE  DISPLAY 


THIS  PROGRAM  DISPLAYS  THE 
DISTRIBUTIONS  ALONG  A  CUT 
AND  THROUGH  THE  ELECTRODE 
GRAPH  ON  THE  CALCOMP  PLOTTER 
SCHOOL  SL6ROUTINE  DRAWP.   IT 
PREDETERMINED  PARAMETERS: 


POTENTIAL  AND  CHARGE  DENSITY 
PERPENDICULAR  TO  THE  WALL 

IT  DISPLAYS  THE  RESULT  AS  A 
USING  A  NAVAL  POSTGRADUATE 
REQUIRES  THE  FOLLOWING 


P  POTENTIAL  DISTRIBUTION 

G  ELECTRON  DENSITY  DISTRIBUTION 

SI  ICN  DENSITY  DISTRIBUTION 

RAT  CHARGE  DENSITY  BOUNDARY  VALUE 

NM  ARRAY  LENGTH 

NC  ELECTRODE  LOCATION 

ZL  STEP  SIZE  MULTIPLE 

TITLE  DESIREO  GRAPH  LEGEND 


OTHER  PARAMETERS  ARE  EXPLAINED  IN  SUBROUTINE  DRAWP 

IMPLICIT  REAL-8  (A-H,0-Z) 

INTEGERS    IT3<12)/12*0/ 

REAL*4    RTB( 28 1/28*0.0/ 

REAL*4    LABL(3)/"P0T    «,«POS    «,'NEG     '/ 

DIMENSION    X(51),Y(51),Z(51),W(51) 

D  I  KENS ICN    A (51, 51) 

EQUIVALENCE     ( TITLE ,RTB ( 5)) 

REAL*8    TITLE( 12) 

READ(5,101)  TITLE 
SET  CCNSTANTS 

RAT=l.D-03 

NM  =  51 

NC  =  26 

ZL=5.D  C2 

H=ZL/CFL0AT<NM-1 ) 

ITB(4)=10 

R  T  B  (  2  )  =  .  1 2  5 

ITB( 5)=1 

ITB(6)=1 

ITB( 8)=2 

DO  12  I=1,NM 
12  X(  I)=H*DFLOAT(  1-1) 
INPUT  DATA 

READ(l)  A 

AA=A(1 ,26) 

DC  1  1=1, NM 

1  YU)  =  A(  I,NC)/AA 
READd  )  A 

DC  2  I=1,NM 

2  Z( I)=Ad,NC) 
READd)  A 

DO  3  1=1, NM 

3  W(I)=A(I,NC) 
NFLAG=1 

10  GO  TO  14,5,6,7) ,NFLAG 

5  DO  8  1=1, NM 

8  Y(I)=MI)/RAT 
GO  TO  4 

6  DC  9  1=1, NM 

9  Y( I)=2(I)/RAT 

4  RT6(3)=LABL(NFLAG) 
ITB( 1 )=MFLAG 

OUTPUT 

CALL  DRAWP(NM,X,Y,ITB,RTB) 

NFLAG=NFLAG+1 

GO  TO  10 

7  STCP 

101  FCRMAT(6A8) 
END 
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c  *  * 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ELECTRON  TEMPFRATURE  PROFILE 


## 


HIS  PROGRAM  CALCULATES 
LONG  A  CUT  PERPENDICUL 
LFCTRCOE.  IT  DISPLAYS 
ALCCMP  PLOTTER  USING  T 
UBRCUTINE  DRAkP.  IT  R 
RE-DETERMINED  PARAMETE 


P  POTENTIAL  DIST 

G  ELECTRON  DENSI 

SI  ION  DENSITY  DISTRIBUTION 

GAM  JOULE  HEATING  PARAMETER 

Z  STEP  SIZE  MULTIPLE 

NM  ARRAY  LENGTH 

NC  ELFCTRCDE  LOCA 

TITLE  DESIRED  LEGEND 


THE  ELECTRON  TEMPERATURE  RATIO 
AR  TO  THE  WALL  AND  THROUGH  THE 

THE  RFSULT  AS  A  GRAPH  ON  THE 
HE  NAVAL  POSTGRADUATE  SCHOOL 
EUUIRES  THE  FOLLOWING 
RS  J 

RI3UTION 
TY  DISTRIBUTION 


1 
C       C 


10 
10 


THER    PARAMETERS    ARE    EX 

IMPLICIT    REAL-8    ( A-H 

INTEGERS    I  TB(  12  J  /l  Z 

REAL-4    RTB( 281/23*0. 

DIMENSION  P( 51  ,51) ,G 

EQUIVALENCE  (TITLE, R 

REAL* 8  TITLEC12) 

READ(5,101)  TITLE 
ET  CONSTANTS 

GAM=1.0D  03 

Z=5.CC  02 

NM  =  51 

NRCM=50 

H=Z/DFLOAT(NROM} 

NC  =  26 

NCM=NC-1 

NCP=NC+1 

G1=GAN*.25D  00/(H*H> 

G2=GAP/(H*H) 

ITB(4)=10 

ITB( 5)=1 

IT6( 6)=1 

ITB( 8)=2 

READ(l)  P 

READil)  G 

PEAD(l)  SI 

DO  12  1=1, NM 
2  X(  I)=H*PFLOAT(  I-i) 

TH(1  )=i  .D  00 

TH(NM)=1.D  00 
ALCULATE  ELECTRON  TEMP 

DO  1  I=2,NRGM 

FL=G(  I,NC) 

DPX=P(I+1  ,N'C)-P(  1-1, 

DPY=P(I,NCP )-P(I ,NCM 

DNX=b< 1  +  1 ,NC) -G< 1-1 , 

DNY=G(I ,NCP)-G(I ,NCM 

ARG=G1*(DNX*DPX+DNY* 

ARG2=G2*( DPX*DPX+DPY 

TH(I )=(-ARG+DS£RTl AR 
1  CONTINUE 
UTPUT  RESULTS 

WRITE (6,  ICO)  (TH(I)f 

CALL  CPAWP(NM,X,TH,I 

STOP 

FORMAT(5H3X,F12.3,/ 

F0RMAT(6A8J 

END 


TION 

CN  GRAPH 
PLAINED  IN  SUBROUTINE  DRAWP 
,0-Z) 
*0/ 
0/ 

(5  1,51),SI(51,51),TH(51),X(51) 
TB(5) ) 


ERATURE  RATIO 


NC) 

) 

NC) 

) 

DPYi/EL-l.D  00 

*DPY) 

G*ARG+ARG2) )/2.D  00 


1  =  1  ,NM) 
TB,RTB) 

)  ) 
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c** 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CONTOUR  PROGRAM 


** 


HIS  PRO 
HARGE  D 
HE  SAME 
OLLCWIN 

AD 

NRCWS 

NXCLS 

ML 

NL1 

IW 

IH 

TABL 
(  1 
(2 
(3 
(4 
(5 
( 

LTG 
( 


GRAM 

ISTR  I 
CAT  A 

G  PAR 
PCTE 
ARK  A 
ARRA 
NUMB 
NUMB 
CESI 
DESI 
6  FO 

) 


5) 
6) 


JOU 
TEM 
ELE 
CHA 
HAL 
CHA 
3  GR 
)  TRU 


(2)  TRUE 


(3 


HE  PROG 
MA 

FL 

CL 

CC 


FAL 
)  TRU 


FAL 
FAM  H 
IN... 

CF... 

VL..  . 

NTUR. 


N/BL6 

N/BL1 

8  TIT 

i 

i 

i 

VCLTAGE  • 


AP   ■ 
S  TIT 
8  TAB 
8  AD( 
AL*1 
SICN 
S  I  ON 
IST/T 
IST/O 
5, TAB 
5,DOL 
(IW*1 
(  IH*1 


10 


CGMMC 
COMMC 
REAL* 

1« 

2' 

3' 

4» 

5' 

6« 

7»  AL  M 
REAL* 
REAL* 
REAL* 
LOGIC 
DIMEN 
DIMEN 
NAMEL 
NAMEL 
REA0( 
READ( 
XVAL  = 
YVAL  = 
N1=0 
DO  10 
Nl  =  Nl 
X  L I N  < 
YLIN( 
DO  1 


PLOTS 

BUT  I  ON 

CELL 
AMETER 
NTIAL 

Y  LENG 

Y  WIDT 
EP  OF 
ER  OF 
RED  OU 
RED  OU 
RMATTE 
LE  HEA 
PEPATU 
CTROCE 
RGE  DE 
L  PARA 
RACTER 
ID  SEL 
E  •  • • AL 

IN 
RE 
SC 
FR 
DM 
E ...  RE 
LI 
GR 
SE.  .CM 
AS  THE 
GENERA 
OUTPUT 
ORIENT 
LABELL 
DETERM 
PLOTTE 
. .  LABE 
ALLOWS 
PLOT. 
SUCH  A 
AMD  IS 
/NROWS 
1/XLIN 
LE(30) 
? 

T 
» 
I 


THE  CONTOURS  OF  THE  POTENTIAL  A 

S  USING  THE  CALCOMP  PLOTTER.   I 

AS  THE  SHEATH  PRUGRAM.   THE 

S  MUST  BE  PREDETERMINED: 

OR  CHARGE  DENSITY  ARRAY 

TH 

H 

POTENTIAL    CONTOURS     TC    BE    PLOTTE 

CHARGE    DENSITY    CONTOURS    TO    BE    P 


IN    INCHES    <    10 
IN     INCHES    <    16 


SE 


TPUT  WIDTH 
TPUT  HEIGHT 
D  VARIABLES 
TING  PARAMETER 
P  P 

"potential 

nsity  boundary  condition 

METER 

1ST  IC  LENGTH 
ECTIONS  TRUE/FALSE 
L  EXTERIOR  CONTOUR  SEGMENTS  AND 
TERNAL  MAXIMUMS  WILL  BE  LABELED 
QUEST  TIC  MARKS  AND  CORRESPONDI 
ALE  VALUES  ONE  INCH  APART  ON  EX 
AME  OF  THE  CONTOUR  GFAPH 
IT  TIC  OPTION 

QUEST  A  ONE  INCH  BY  ONE  INCH  ST 
^,E  GRID  SUPERIMPOSED  ON  THE  CON 
APH 

IT  GRID 
FOLLOWING  SUBROUTINES: 
L  COORDINATION  AND  PREPARATION 


ND 

T  USES 


D  <  21 
LOTTED 


S  ARR 
ING  C 
INES 
D 

LLING 

SUPE 

IT  U 

S  RES 

A  MO 

tNCOL 

(100) 

/« 


AYS  TO  CORRESPOND  WITH  CUT 

ONVENTION 

LINES  CF  CONSTANT  VALUE  TC 

AND  PLOTTING  COORDINATION 
RPOSITICN  OF  3  CONTOURS  ON 
TI LIZES  ADDITIONAL  SUBROUT 
TOFt  SCAM,  TRACE,  CALC  AND 
DIFICATION  OF  AN  NPS  PRCGP 
S  f  I  W  » I H ,  N 1 
,YLIN( 100) 


NG 
TERIOR 


RAIGHT 
TOUR 


CF 
PUT 
BE 


ONE 
INES 
PLOTT 

AM. 


TEMP 


R.  DOLS 


LEA(4) 

L  E  (  6  ) 

51,51) 

LTG(3) 

A(81  ,6 

CL(30) 

ABL/TA 

OLL/NR 

L) 

L  ) 

.  )  /  {  NC 

.)/(NR 


•DENSITY 

i 

i 

'LENGTH 
•ON 

•BOX  1219 
/•  ELECTRO 


GAM 


•BET 


POTENTI 


'  t' 


ION  CHAR 


•  ,»GE 


1) 


BLEtLTG 

0WS,NC0LS,NL,NL1,IW, IH 


0LS*1.-1 
OWS-l  .-1 


1  =  1,1 

+  1 

M)=XVAL*(20) 

N1)=0.0 

K=l,3 
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19 


11 


3 
12 


INPUT  ARRAY 
REAQ(2)  AO 
DO  19  I=l,NROWS 
DO  19  J=1,NCCLS 
A( I, J)=AD( I  ,  J) 
CALL  FLOP (A  , NROWS, NCOL S ) 
FIRST  FIGURE  CONTOUR  LEVELS 
CALL  CLVL( A,CL,NL) 
WRITE(6,210) 
DO  3  1=1, NL 
WRITE(6,2i5)  I,CL(I) 
TI 7LE(8}=TABLE(1) 
TITLE( 11) =TA6LE(4) 
TITLE( 14)=TA6LE( 2) 
TITLE( 17)=TA8LE(3) 
TITLE(20)=TA3LE( 5) 
TITLE( 23)=TABLE(6) 

CALL  CCNTUR ( A , N'ROWS , NCOLS , 8 1 ,CL ,NL 
MOW  FLCT  CHARGE  FIELD  CONTOURS 
I-F(K.EQ.31  STOP 
KL=2*K-1 
KF=KL+1 

TITLE(27) =T  ITLEA(KL) 
TITLE(28)=TITLEMKP) 
CONTINUE 

FORMAT  (•     ' ,T26,'THE 
FCPMATC     '  ,  T3  3,  'CL(  » 
STOP 
ENC 


1 

210 
215 


TO  BE  PLOTTED  FOR  MATRIX  A 


TITLE, LTG) 


LEVELS  PLOTTED' ) 
,12, l  )=»  ,F8.4) 


SUBROUTINE  FLOP( Z , MRCWS , NCOLS J 


D  I  FENS  ION  Z(81 
IIVRT=NRuWS/2 
DO  3  1=1,1 IVRT 
M=NRCkS-{I-l) 
DO  3  J=l, NCOLS 
SAVE=Z( I, J) 
Z( I, J)=Z(M, J  ) 
Z(M,  J)  =  SAVE 
RETURN 
END 


61) 


CM  I N= CM (I 
CMAX=CM(I 


SUBROUTINE    CLVL( CM, C LM , NU ML ) 

C0MMCN7BL6/NR0WS ,  NCOLS , Id , IH,N1 

DIMENSION    CM( 81,61) , CLMCNUML) 

CMIN=CM< 1,1 ) 

CMAX=CMIN 

DO    5    J=l, NCOLS 

DO    5     I=1,NkCWS 

IF (CM* I, J) . LT.CMIN) 

IF(CM< I, Ji .GT.CMAX) 

CONTINUE 

NOk    DETERMINE    CONTOUR 

J=NUML 

I=NUML-1 

ANL=I*1. 

PLIMT=(CMAX-CMIN)/ANL 

NOW    FILL    THE    CONTOUR    LEVEL 

DQ    6     1=1, J 

CLM  I)=CMIN*(  1-1  )*PLINT 

RETURN 

END 


LEVELS    TO    BE    PLOTTED. 


VECTOR 


SUBROUTINE  C ONTURC AM ,M , N, MX ,CL , NL, T I TLE , LTG ) 

COMMON  /INTFAC/  X,Y 

COMMON  /DAYFOF/  MT , NT , N  I  ,  IX ,  I Y  ,  I  DX ,  I DY,  I SS ,  I T,  I  V, NP , 
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1 

60 
71 

64 


2 

40 
61 


2011 


C 
C 


CO 

CO 

CO 

CO 

CO 

RE 

RE 

DI 

DI 

DI 

DI 

LC 

JC 

LA 

CH 

WH 

IF 

WR 

FO 

WR 

FO 

RE 

CH 

IF 

WR 

FO 

IW 

NO 

IF 

WH 

GO 

DI 

DI 

XN' 

YM 

SL 

SL 

DI 

DI 

DI 

DI 

DI 

DI 

DI 

DI 

DI 

DI 

DC 

DI 

01 

ST 

CA 

CA 

CA 

DI 

DI 

DI 

DI 

DI 

DI 

DI 

DI 

DI 

DI 

CA 

SL 

SL 

IE 

IE 

IF 

DR 

ST 


MM  ON 
MMCN 
MMCN 

MMCN 

AL*8 

AL-8 

MENS 

HENS 

MENS 

NENS 

G1CA 

=  0 

BL  =  L 

ECK 

ICH  = 

(  IW) 

IT  E  ( 

RMAT 

1TE( 

hMAT 

TURN 

ECK 

(  IW- 

ITE( 

RMAT 

=  9 

W  CH 

(  1H) 

ICH  = 

TO 
TSDX 
TSOY 
IN=1 
IN=- 
CPEX 
CPEY 
TSX( 
TSX( 
TSX( 
TSX( 
TSX( 
TSY( 
TSY( 
TSY( 
TSY( 
TSY( 

201 
TSX( 
TSY( 
ARTP 
LL  F 
LL  P 
LL  L 
TSX( 
TSX( 
TSX( 
TSX( 
TSX( 
TSY( 
TSY( 
TSY( 
TSY( 
TSY( 
LL  L 
GPEX 
CPEY 
NDX  = 
NDY  = 
(  .NO 
AW  T 
ART 


/TA3L/    TA6C(20,6) ,  JC 

/DITS/XMIN,YMIM,SLOPEX,SLOPEY,DITSDX,DITSDY, 

/BL6/i\R0v,S,NCDLS,  IW,  IH,N1 

/3L 11 /XL  IN (100)  ,YLIN( 100) 

/BL13/NSCN 

TITLEQJ 

WIDTH/ 'WIDTH* / , He IGHT/ • HE IGHT* / , WHICH 
ICN  AM(MXtl) ,CL(  1) 

ION  REC(900),   X<1800),  YU800) 

ION  IPT(3,3)  , INX( 8)  ,  INY(8) 
ICN  DITSX(5) ,DITSY( 5) 
L*l  LTG( 1) , MINUS, LABL 

TG(1  ) 

IW    PARAMETER 
kIDTH 
1»1,  2 

6,60)     WHICH 

CO' ,T7,A8,'0F    CONTOUR    GRAPH     ILLEGAL.') 


«0' 


T7,'N0  GRAPH  WILL  BE  PRODUCED.') 


IF  IW  IS  TOO  WIDE 

9)  3,3,40 

6,61) 

('0«,T7,'IW  PARAMETER  GREATER  THAN  9.  CONTUR 


0/ 
0/ 


ECK  I 
4,4, 
HEIGH 
1 

=  <N-1 
=  (-1  . 
.0 
M 

=  1 
=  1 

1)  =  1. 

4)  =  1. 
5)=1. 

2)  =  N 
3)=N 
1)=-1 
2)=-l 
5)=-l 
3)=-M 
4)=-M 
1  1=1 
I  J=SL 
I  )  =  SL 
=  (9.0 
LOTS 
LOT(S 
INE1D 

1)  =  DI 
5)=DI 
4)=DI 

2)  =  DI 
3)=DI 
1)=DI 

5)  =  DI 
2)=DI 

3)  =  DI 
4)=DI 
INE(D 
=  1.0/ 
=  1.0/ 
SLOPE 
SLOPE 
T.LTG 
IC  MA 
CN  LF 


PARAMETER 


.0)/IW 
0+M)/IH 


DITSDX 
CITSCY 

0 


.0 
.0 

.0 


,5 

OPEX*(DITSX( I )-XMIN) 
OPEY*(DITSY( I )-YMIN) 
-  I W  )  /  2  . 0 

TARTP,0.0 ,-3) 

ITSX,DITSY,5,1,1) 

TSXU  J-.5 

TSX(1 ) 

TSX(4)-.5 

TSX(2)+.5 

TSX(3 )+.5 

TSY(1  )+.5 

TSY(  1  ) 

TSY(? )+.5 

TSY(3)-.5 

TSY(4J-.5 

ITSX,DITSY,5,1,1) 

DITSDX 

CITSDY 

X*N  +  1 

Y*M  +  1 

(2)  )  GO  TO  34 

RKS  ON  OUTER  FRAME 

FT  EDGE  GOING  DOWNWARD 
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2022 


21 


,22,23,24) ,  IFLAG 


RIGHT  EDGE  GOING  DOWNWARD 


BOTTOM    EDGE 


C 
C 
C 


IFLAG=0 
ZINGX=-.l 
ZINGY=C.O 
ZX=0.0 
ZY=-.9 
CX  =  DITSX(1  ) 
CY  =  [)ITSY<  1  )-.5 
I  END =11 
2222     IFLAG=IFLAG+1 

DC    2022    1=1,1  END 

CALL    PLCT(CX,CY,3) 

CGCRDX=CX+Z INGX 

CCCRDY=CY+ZINGY 

CALL     PL0T(CG0RDX,C00RDY,2) 

CX=CX+ZX 

CY=CY+ZY 

GO    TO    (21 

NOW    DO    THE 

ZINGX=.l 

CX=DITSX<2) 

CY  =  DITSY( 2)-. 5 

GO    TO    2222 

NCW    DC    TOP    EDGE 

22  ZINGX=C.O 
ZINGY=.l 
ZX=.9 
ZY=O.C 

CX=DITSX( 1  )+.5 
CY=DITSY(1 ) 
I  END* 11 
GO    TC    2222 
NOW    DC    THE 

23  ZINGY=-.l 
CX=DITSX(4)  +.5 
CY=DITSY(4) 
ZINGY=-.l 
GC    TC    2222 
NGU    LABEL    TIC    MARKS 
DO    X-DIRECTION    FIRST,     TOP    EDGI 
POSIT ICN    PEN 
DELTAX=.l 
IFLAG=C 
ZX=.9 
ZY=O.C 

CX  =  0I1SX(  1 J+.35 
CY=DITSY( 1I+.12 
IFLAG=IFLAG+1 
XZERC=0.0 
DO  3333  1=1,11 
CALL  NUM3ER(CX,CY, 
CX=CX+ZX 
CY=CY*ZY 

XZERO=XZERO+DELTAX 
GO  TO  (31 ,32,33,34) , 
LA3EL  BOTTOM  EDGE  TI 
CX  =  DITSX(4H-.35 
CY  =  D  ITSY(4)-.19 
GC  TO  3033 
LABEL  LEFT  EDGE 
CX=DITSX(4)-.4 
CY=DITSY<4)+.46 
DELTAX=.l 
IEND=IENDY 
ZX^O.O 
ZY  =  .9 

GO  TO  3033 
NOW  LABEL  RIGHT 

33  CX=D  ITSX( 3 )  +  .12 
CY=DITSY(3)+.46 
GC  TO  3033 
CHECK  IF  GRID  1ESIRFD 

34  CALL  PESTOF(LTG, IENDX, IENDY,NL, AM, M, N,MX , CL , STARTP,TIT 


24 


3033 


^Jjj 


31 


32 


14,XZERO,0.0,2) 


IFLAG 

MARKS 


OF  TIC  MARKS 


EDGE  TIC  MARKS 
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4044 


4444 


4  2 


35 


20 


777 
778 


779 


LE,DITSX,    DITSY) 
RETURN 
END 
SUBROUTINE  RES TOF( LTG , I ENDX , 

TPtTITLEt    DITSX, DITSY) 
REAL*3  TITLE(l) 
DIMENSION  AM(MX,1)  ,CL<  1  ) 
DIMENSION  OITSX( 5) ,DITSY( 5) 
L0GICAL*1  LTG( I) , MINUS, LABL 
C0«MCN/TABL/TA3C(20,6) ,  JC 
CO WMON/D I T  S/X  M I N , YM I N , S  LQ  PEX 
C0MM0N/BL6/NR0WS,NCGLS,IW,IH,N1 
COMMCN/BLll/XLINdOO)  ,  VLIN{  100) 
C0MMCN/BL13/NSCN 
IF( .NCT.LTG(3I  )  GO  TO  35 
DRAW  INCH  BY  INCH  GRID 
IEl\D=IBNDX-2 
POSITION  PEN 
IFLAG=C 

CX=DITSX(1 J +.5 
CY  =  DITSY(1  )-.5 
COCRDX=0.0 
CCCRDY=-IH 
DX  =  1  .0 
DY=0.0 

DO  4444  1=1 , IEND 
CX=CX+DX 
C  Y=C  Y  +  CY 

CALL  PL0T(CX,CY,3) 
ZX=CX+CGCRDX 
ZY=CY+CCORDY 
CALL  PLCT(ZX,ZY,2) 
IF(IFLAG)  35,42,35 
IFLAG=1 
IE.\D=IENDY-2 
CY  =  D  ITSY(4)  +  .5 
CX=DITSX(4)+.5 
CCCRDX=IW 
CCCRDY=0.0 
DX=0.0 
DY  =  1  .0 
GO  TO  4044 
CONTINUE 
CALL  L  I N  E  (  X  L  I  N ,  Y  L  I N  ,  N 1  , 1 ,  -  6  ) 


IENDY,NL,AM,M, N,MX,CL,STAR 


SLOPEY,DITSDX,DITSDY, 


M ,  N  ,  MX  ,  C  L  (  I )  ) 


DO  2  0  1=1 , NL 

CALL  SCAN(AM, 

NSCN=1 

IF(.NOT.LAbL)  GO  TO  778 

IF(JC.EQ.O)  GO  TO  778 

DO  777  1=1, JC 

COCRDX=TABC< I ,4) 

COCRDY  =  TABC(  1,5) 

CLEV=TA8C( I ,6) 

CALL  NUMB ER ( COGR DX , COQRDY  , 

CONTINUE 

CALL  SYMBOL (-STARTP 

CALL  SYMBOL (-STARTP 

CALL  SYMBOL (-STARTP 

CALL  SYMBOL (-STARTP 

CALL  SYMBOL (-STARTP 

CALL  PLCT(-STARTP,IH+6.5,-3) 

CALL  PLOTE 

RETURN 

CALL  SYMBOL(-STARTP, IH  +  1.0,  .21 

CALL  PLOT(-STARTP, IH+4.5,-3) 

CALL  PLCTE 

PETUPN 

END 

SUBROUTINE  SC  AN'(  AM,  M  ,  N  ,MX  ,CL  ) 

DIMENSION  AMtMXtl  ),REC(900), 

DIMENSION  IPT(3, 3) ,INX(8) 

COMMON  /DAYHOF/  MT, NT, NI, 


07,CLEV,O.O,3) 


IHH.O, 
IH+1.5, 
IH+2.0, 
IH+2.5, 
IH+3.0, 


.21 
.21 
.21 
.21 
.21 


TITLE(25) ,0.0,48) 
TITLE  (19) ,0.0,48) 
TITLE(  13)  ,0.0,48) 
TITLE (7) .0.0,48) 
TITLE  (1  )  tO.Ot 48) 


TITLE(25) ,0.0,48) 


X(1300),     Y(1800) 
, INY(8) 
IX, I Y, IDX, IDY, ISS, IT, I V,NP,NC 
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53 


55 
57 


110 


15 
17 


»JT, 

IKY,DL,RA,THE 
iCN    /INTFAC/    X,Y 


PY,REC,CV,     IFT,INX, 


COMMON  /INTFAC/  X,Y 
LOGICAL  *1  LABL, MINUS 

COMMON/DITS/XMIN,YMIN,SLOPEX, SLOPEY,DITSDX,DITSCY, 
C0MM0N/BL13/NSCN 
D  =  0. 


R=l. 

TH  = 

NP=0 

DL=D 

RA=R 

TFE=TH 

MT=N 

NT=M 

CV=CL 

IF(NSCN.EQ.l)  IZW=0 

I F(I ZW-120631)  1,3,1 

I PT(  I*  1)=8 

I  P  T  (  1 ,  2  )  =  1 

!OT!  I  .  -j  \  -  -> 


x   r  I  \  L  1  C  I  —  L 

IPT(  1  ,3)=2 
IPT(2,1}=7 
IPT( 2,3)=3 

TOT/3. 1\-A 


IPT(3,1)=6 

IPT( 3,2)=5 

IPT(3,3)=4 

INX(  1)=-1 

INX( 2)=-l 

INX( 31=0 

I  N  X  (  4  )  =  1 

INX( 5)=1 

INX(6)=1 

INX( 7)=0 

INX( 8)=-l 

INY( 1 )=0 

INY( 2)=1 

INY(3)=+1 

INY(4)=+1 

INY( 5)=0 

INY( 6)=-l 

INY( 7)=-l 

INY< 8)=-l 

IZW=120631 

XT=MT 

DO  58  J=l,900 

REC( J)=0 

ISS=0 

MT1=MT-1 

IDIP.=  1 

DO  110  1=1, MT1 

IF(AM(  1,1  )-CV) 

IF(AM<1,  1+1  )-CV 

IX=I*1 

IY  =  1 

IDX=-1 

IDY=0 

CALL  TRACE  (AM, MX) 

CONTINUE 

NT1=NT-1 

IDIR=2 

DO  20  1=1, NT1 

IF(Atf<  I,HT)-CV)  15,2  0,20 

IF(AMI+1,MT)-CV)  20,17,17 

IX  =  MT 


55,110, 110 
>  110,57,57 


IX=MT 
IY=I+1 

IDX  =  0 
IOY=-l 
CALL  TRACE 
20  CONTINUE 
r  n  t  p=-a 

22 


(AM, MX) 


IDIR=3 

DO  3  0  I=1,MT1 
MT2=MT+1-I 
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25 
11 


30 


35 
37 


40 


5 
7 

12 

9 
11 


10 


17 

18 
5 


47 


49 


IF(A 
IF(A 
IX  =  M 
IY  =  N 
IDX  = 
IDV  = 
CALL 
CGNT 
IDIR 
DO  4 
NT2  = 
IF(A 
IF(A 
IX  =  1 
IY  =  N 
IDX  = 
IDY= 
CALL 
COM 
IDIR 
ISS  = 
NT1  = 
KT1  = 
00  1 
DC  1 
IF(A 
IF(A 
CGK  = 
IF  ( 
00  9 
IF  < 
CGNT 
IX  = 
IY  =  J 
IDX  = 
IDY  = 
CALL 
CO  NT 
RETU 
END 
SUBR 

0  I  KE 

DII^L 

CO  KM 
CO  KM 
PY  =  0 
RC  = 
RS  = 
JT  =  0 
N  =  0 
IXC  = 
IYO  = 
ISX= 
ISY  = 
IS=I 
JTB  = 

1  S0  = 
IF(  I 
ISC  = 
IT=0 
CO  NT 
CALL 
NZ  =  N 
N  =  NZ 
IF    ( 

xs  =  x 

YS  =  Y 
X(N- 
Y(N- 
X(N) 
Y(N) 
15  =  1 


K( 

K< 
T2 
T 
1 

0 

T 
IN 
=  4 
0 

NT 
M( 
Ml 


NT,MT2)-CV)     25,30,30 
NT,MT2-1J-CV)     ^0,27,27 


RACE 


(AM, MX) 


1=1 , NT1 

+  1-1 

NT2,1)-CV)    35,40,40 

NT2-l,i)-CV)    40,37,37 


T2- 

0 
1 

TP 
IMJ 
=  5 
1 

NT- 
MT- 
0    J 

0  I 
M(  J 
M(  J 
ICO 
NP) 

ID 
REC 
INL 

1  +  1 


I 


ACE     (AM, MX) 
E 


2,NT1 

1,MT1 
I  )-CV)5, 
1  +  1  )-CV) 
(  1  +  1  )+J 
12,  11,  12 

1,  NP 
IDJ-COM) 


10,10 
10,7 


9,10,9 


-1 
0 

TRACE 

INUE 
RN 


(AM, MX) 


OUTIME  TRACE  (AM, MY) 

NSIGN  AM(MY,1  ), REC (900),   X(i800), 

NSICN  IPT(3,3),  I N  X ( S ), INY(8) 

C-N  /OAYHOF/  MT,NT,NI,  IX,  IY,  IDX,  IDY, 

ON  /INTFAC/  X,Y 

.0 

CCS  (THE)^RA 

SIN  (THE)*RA 


IX 
IY 

IOX  +  2 
IDY  +  2 
PT( ISX, 
0 

IS 

SC-8)  18 
I  SC-8 


Y( 1800) 
ISS, IT, IV,NP, 


ISY) 


18,17 


INUE 
CALC 


IT+JT-1 

(i^-n 

(N-l) 
1 )=X(M) 
1)=Y(N) 

=  xs 

=  YS 
S  +  l 


(  A  M  ,  MY  ) 


)  49,49,47 
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308 
103 

51 

20 
21 
22 
23 


10 
13 
19 
11 
12 
206 
213 

217 

214 


6 

306 

16 


50 

307 

73 

74 


41 


JT=I 
IF  ( 

IS=I 
IDX  = 
IDY= 
1X2= 
IY2  = 

jTe= 

IF  ( 

PR  IN 

FORM 

RETU 

CCNT 

IF  ( 

I  F  (  I 

I  F  (  I 

I  F  C  I 

CONT 

CALL 

GO  T 

IF(  I 

IF  ( 

IF  { 

IF  ( 

IF(C 

IF  ( 

DCP^ 

IF  ( 

IF  ( 

I  X  =  I 

IDX  = 

PY  =  2 

CALL 

IX=I 

GC  T 

IY  =  I 

IDY= 

PY  =  2 

CALL 

IY  =  I 

IF(A 

KP-U 

REC( 

IS  =  I 

IX=I 

IYM 

GO  T 

XT=M 

IF(A 

NP=N 

REC( 

DO  7 

X(  I  ) 

Y(  I  ) 

CALL 

RETU 

END 

SUBR 

DIME 

DIME 

coy.M 

CO  MM 
IT=0 
N  =  N«- 
IF  ( 
IF  ( 
X(N) 
Z=IY 
IY2  = 
DY  =  I 
Y(N) 
RETU 


T 

IS-9) 

S-8 

INX(  I 

I  NY  (  I 

I  X  +  I  D 

1Y  +  I0 

JT6  +  1 

JTB-1 

T 

AT(1H 

PN 

INUE 

ISS) 

X-IXO 

Y-IYO 

S-ISO 

INCF 

CALC 
0  73 
X2)  1 
IX2-M 
IY2) 

IY2-N 

V-AM{ 

IDX** 

{  AM(  I 
DCP-C 
I N  X  (  I 
X+IDX 
-  IDX 
.0 

CALC 
X  +  IDX 
0  6 
Y+IDY 
-IDY 
.  C 

CALC 
Y+IDY 
M.(  IY, 
F+l 
NP)  =  1 
S  +  5 
X2 
Y2 
0  9 
T 

M  (  I Y , 
P+l 
N  P  )  =  1 
4  1=1 
=  X(I  ) 
=  RS*Y 

PLOT 
PN 


8,7,7 

S) 
S) 

X 
Y 

799)  51,51,308 

103,CV,X(N) ,Y(N) 
0,23HA  CONTOUR  LINE  AT  LEVEL,E12.5, 


10,10,20 
)  12,21,12 
)  12,22,12 
)  12,23,12 

(AM, MY ) 

3,50,13 

T)  19,19,50 

11 ,50,11 

T)  12,12,50 

IY2,IX2))  206,206,5 

2+ ID Y** 2-1)  213,6,213 

Y,  IX)+AM(  IY,  IX2)+AM(  IY2,  IX)  <-AM(  IY2,  1X2  )  )/4.0 

V)  5,217,217 

S-l))  214,215,214 


(AM, MY) 


(AM, MY) 
IX-D-CV)  306,16,16 
00*IX+IY 


IX-D-CV)  307,73,73 

CO*IX+I Y 
,N 

+RC*Y(I ) 
(  I  J 
T(N,CV) 


OUTINE  CALC(AM,MY) 

NSION  AM(MY,1  ),REC(900),   X(1300),  YU800) 

NSION  IPT(3,3J ,INX(8) ,INY(8) 

ON  /DAYHOF/  MT , NT , N I ,  IX , I Y ,  I  DX ,  I  DY , I SS , I T , I  V , NP , N, 

ON  /INTFAC/  X,Y 


] 

ICX**2  +  IDY**2 
ICX)  10,2,10 
=  IX 


-1)  20,1,20 


IY+IDY 

=((AM(IY, IX)-CV)/( AM( IY,IX)-AM( IY2,IX) ) )*DY+Z 
RN 
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10 

44 
20 


24 
21 
23 
27 


25 

33 

28 


100 


200 


Y(N) 

W=IX 

DX  =  I 

1X2  = 

X(N) 

RETU 

1X2= 

IY2= 

W=IX 

Z=IY 

DX  =  I 

DY  =  I 

DCP  = 

IF  ( 

IF  ( 

AL  =  A 

V=.5 

X(N) 

Y(N) 

PY  =  0 

RETU 

IT  =  1 

AL  =  A 

V=.5 

X(N) 

Y(N) 

Y  (  N ) 

PETU 

ENC 

SUBR 

COMM 

LCGI 

CO  MM 

CO  MM 

SCAL 

DO  1 

X(  I  ) 

Y(  i  j 

CALL 
IF(. 
DIR= 
GO  T 
DIR  = 
CCCR 
COCR 
CALL 
RETU 
MOVE 
DIR  = 
CCCR 
CCCR 
GC  T 
MOVE 
CCCR 
COCR 
GO  T 
SEAR 
XMAX 
SMIN 
YMIM 
YMAX 
VM  IN 
DO  2 
IF(X 
IF(Y 
IF(Y 
IF(X 
SMIN 
YMIN 
CO  NT 
JC  =  N 
IF(J 


=  IY 

DX 

IX+IDX 

=((AM( IY, 

RN 

IX+IDX 

IY+IDY 


IX)-CV)/< AM( IY, IX)-AM( IY, 1X2) ) ) *DX*W 


DX 

DY 

(  AM(  I 

PY-2* 

DCP-C 

M  (  I  Y  , 

*<  AL  + 

=  V*DX 

=  V*DY 

.0 

RN 


Y, IX)+AM( IY 
0)  24,21,24 
V )  21,21,25 
IXJ-DCP 
DCP  -CV)/Al 
+  W 
♦  Z 


1X2) +AM(  IY2, IX)  +  AM(  IY2, 1X2  )) /4.0 


M( IY2, 1X2) -DCP 
*(AL+DCP-CV)/AL 
=-V*DX+W  +  DX 
=-V*DY+Z  +  CY 
=-V*DY«-Z  +  DY 
RN 


CU 
ON 
CA 
ON 

ON 

E 

OC 

=  s 

=  s 

L 
NC 
0. 
0 

90 
DX 
DY 

N 
RN 

P 
90 
DX 
DY 
0 

P 
DX 
DY 
0 

CH 
=  X 
=  X 
X  = 
=  Y 
=  Y 
OG 
(  I 
(  I 
(  I 
1  I 
=  X 

x  = 

IN 
UM 
C) 


TINE 
/INTF 
L*l  M 
/TA3L 
/DITS 
POINT 
1  =  1, 
LOPFX 
LOPFY 
INE(X 
T.LAb 
C 
(1,2, 


PLOTT(NP,CV) 

AC/XI  1800)  , Y(  1800) 

INUS,LABL 

/  TA3C(20,6) ,JC 

/XMIN,YMIN,SLOPEX,SLCPEY, DITSDX,DITSDY, 

S    FOR    PLOT    ROUTINE 

NP 

*<X(I J-XMIN) 

*(-Y( I) -YMIN) 

,Y,NP,1,1 ) 

L)     RETURN 

3,4,6),     IDIR 


=  X(1     ) 

=  Y(  1     ) 

UMBER (COQRDX,COORDY  ,.07,CV,DIR,3) 

EN  DCWN  ONE  HALF  INCH 

=  X(1  ) 

=Y(1  )-.3 

5 

EN  TO  THE  LEFT 

=X<1  )-.3 

=  Y(1  ) 

5 

FOR  XMAX, XMIN,YMAX, YMIN, AND  SAVE  YMINX 
(1) 
MAX 

Yd ) 

MINX 
MINX 

1=2, NP 

.GT.XMAX)  XMAX=X(I) 

.LT.VMIN)  VMIN=Y(I) 

.GT.YMAX)  YMAX=Y(I) 

.GE.SMIN)  GO  TO  200 
(I) 
Y(  I  ) 
UE 
BER  0 

40C, 


F  ENTR IES  IN  TABC 
500,400 
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CONTOUR    TO    BE     INTERIOR    TO    ANOTHER 


1 )  .AND.YMAX.GT.TABC(  1,2) .AND.VMIN. 


400  DO  900  1=1, JC 

IF(XMAX.LT.TABC(  I ,1 )  . AND  .  YMAX . L T . T A BC < I ,  2)  .  AND. VM IN. 
900  CONTINUE 
C      DID  NOT  FIND  THIS 
C      CHECK  IF  EXTERIOR 
DO  1CC0  I =1 , JC 
IF(XMAX.GT.TA6C< I 
1000  CONTINUE 
500  IF  (JC.EQ.20)  RETURN 
JC=JC+1 
MC  =  JC 
600  TABC(MC, 1)=XMAX 
TAEC(MC,2 )=YKAX 
TAEC( KC,3)=VKIN 
TABC(MC,4)=SMIN 
TA8C<MC,5)=YMINX 
TABC(MCt6)=CV 
RETURN 
C      CHECK  IF  THIS  INTERIOR  ONE  IS 
700  IF(CV.LE.TABC(I,6)  )  RETURN 
2000  HC=I 

GO  TC  600 
C      CHECK  IF  LEVEL  OF  THIS  EXTERIOR 
800  IFCCV.LT.TABC.  (  1,6)  )  RETURN 
GO  TO  2000 
END 


OF  HIGHER  LEVEL 


ONE  IS  HIGHER 
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3.   Boundary  Layer  Program 

The  program  which  solves  equation  A.  18  from   Appendix   A 

is  shown  in  the  following  pages.   It  calculates  the  solution 

to  the  boundary  layer  voltage  drop  for   a   varyincr   C  .    It 

P 

therefore  requires  an  enthalpy  table  for  the  calculation, 
and  such  a  table  is  provided  in  data  format  following  the 
program.  This  particular  table  is  for  a 
Toluene/oxygen/cesiura  mixture  at  temperatures  from  1500°K  to 
26U0°K.  An  explanation  on  the  use  of  the  program  is 
included. 
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C   THIS  ROUTINE  I 
C   LAYER  VOLTAGE 
C   SIMPSON'S  RULE 
C   ENTHALPY  TABLE 
C   THE  FCLLCWING 
C      POT=SFED  10 
C      TFREE=FPEE 
C      TWALL=WALL 
C      TTCP=MAX  TE 
C         GREATER 
C      TBCT=LCWEST 
C         EQUAL  TO 
C      NPTS= NUMBER 
C      MAX=NUfbER 
C   ENTHALPY  TABLE 
C0MMCN/BL1/ 
NTH (100 
READ(5t101) 
WRITE (6, 100 


READ(5,102)  TT 


BOUNDAR 

S  DES  IGNE 

DROP  IN  A 

INTEGRAT 

FOR  THE 

PARAMETER 

NIZAT ICN 

STREAK  TE 

TEMPERATU 

MPERATURE 

THAN  OR  E 

TEMPERAT 

OR  SMALL 

OF  POINT 

CF  POINTS 

MUST  BE 
OELtTFREE 
) 

POTtTFPE 
)  POT,TFR 


Y  LAYF 
D  TO  C 

N  MHO 

ION  TE 

RANGE 

S  ARE 

POTENT 

MPERAT 

RE 

IN  FN 
GUAL  T 
URE  IN 
ER  THA 
S  USED 

I  N  EN 
EQUALL 
tTWALL 


R  ROU 
ALC'JL 
CHANN 
CHNIO 
OF  TE 
REQUI 
IAL  I 
URE 

THALP 

0  TFR 
EN  MI 
N  TWA 
IN  I 
THALP 
Y  INC 
,TBOT 


TINE 

ATE  THE  BOUNDARY 

EL.   IT  USES  A 

UE  AND  REQUIRES 

MPFRATURES  CONS  I 

RED: 

N  EVf  S 


AN 
CERED, 


Y  TABLE.  TTOP  MU 
EE. 

ALPY  TABLE.   MUS 
LL. 
NTEGRATION  ROUTI 

Y  TABLE 
REMENTED 

, ALPH tENTZt TTOP t MAX tNtE 


ST  BE 
T  BE 
NE 


READ( 5,103) 
WRITE (6 f 105 
MAX1=MAX-1 
NN=N-1 
NM=N-2 
ALPH=PCT-5. 
DEL=(TTOP-T 
THfeALL=TdAL 
DELM=(TNALL 
M=DELM 
ENTZ=ENTH(M 
TBCT)/D 
DELM=(TFREE 
M=CELN 
IFtM.LT.MAX 
ENTF=ENTH(M 
GO  TO  5 

4  ENTF=ENTri(M 

TBCTJ/D 

5  ENTHM1  =  ENT 
FAC=ENTHM1* 
C0EF2=FLCAT 
C=FCTN(1.0) 
DELTA=(1.0- 

SIMPSCN'S    RULE 
SUM1=0.0 
SUM2  =  C0 
NPT1=NPTS-1 
NFT2=NFTS-2 

COMPUTATION  FO 
SUM=C 

COMPUTATION  FO 
THET A=THWAL 
DC  10  J=2»N 
SUM1  =  SUKH-F 
THETA=THETA 
CONTINUE 

COMPUTATICN  FO 
THETA=THWAL 
DO  20  J  =  3,N 
SUN2  =  SUM2«-F 
THETA=THETA 
CONTINUE 
VOLT=DELTA* 
VGLT=CCEF2* 
WRITE(6,104 
STOP 

FCFMAT(3F9. 

FCRMAT(2E9. 

F0RMAT(3Xt' 

TEMP=« 


E  ,T*'ALL»NPTStN 
EE,TWALL,NPTS,N 
T,MAX 
(EMTH( J) ,  J  =  l ,MAX) 


10 


20 


101 
102 
100 


8C54E 

B0T)/FLCAT(MAX-1 ) 
L/TFPEE 
-TBOT)/DEL 

+  1  )<  (  ENTH(M  +  2)-ENTH(M  +  l  )  >  *<TW  ALL-  M*D  EL- 
EL 
-TBOTJ/DEL 

I)  GO  TO  4 

AX) 

+1 H  ( ENTH(M+2)-ENTH(M+l ) ) *(TFREE-M*DEL- 

EL 

-ENTZ 

*N 

(N)*TFREE/FAC 

THWALL)/FL0AT(NPTS-1  ) 
INTEGRATION  PATTERN  l,4,2,4tl 


R  LAST  POINT.  FIRST  POINT  IS  ZERO. 

R  EVEN  POINTS. . .THOSE  MULTIPLIED  BY  4 

L+DELTA 

PT1T2 

CTN(THETA) 
+2. -DELTA 

R  ODD  POINTS. ..THOSE  MULTIPLIED  BY  2 

L+2.*DELTA 

PT2»2 

CTN(THETA) 

+2.*DELTA 

<SUM+4.*SUMl+2.*SUM2)/3. 

VOLT-1.0 

)  TWALLTTFREE,THWALL,VOLT 

Ct 219) 

IONIZATION  P0TENTIAL=,tF7.3,3X,,FRFE  STREAM 
,F6.1 ,3Xi  'WALL  TEMP=« t F6 . 1 , 3X , • NP TS=  ■  , 14 »  ■ 
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c 
c 
c 


c 

3. 
26 


103    FCPKAT(7E9.0) 

10  4    F0RMAT(3X.2F9.0,F9.3,E12.4) 

10  5    FORMAT (/,7X, • TWALL'  ,4X,'TFRcE» ,  3X, ■THWALL*  ,  3X,  • VOLT*  ) 

END 

FUNCTIGN    THLTA 

FUNCTION    FCTN(THETA) 

CCPMCNVRLl/  DEL,  T  FREE,  TWALL,  T60T , AL PH , ENTZ , TTOP  ,  ,MAX,N,E 

NTH (100) 
MAX1=MAX-1 
MAXM=l«AX-3 
NN=N-1 

A=THETA**(-1.75J 
B=EXP(ALPH*( 1.0/THETA-l.O) ) 
TEPP=THETA*TFREE 
DELMM TEMP-T30T)/DEL 
M=DEL^ 

IF(M.LT.MAX1 )     GC    TO    4 
ENT=ENTH(MAX) 
GO    TO    5 

4  ENT=Ef\TfU M  +  i )+(ENTH( M*2 ) - ENTH( M+l ) ) *( TEMP-M*DEL- 

TBCD/DEL 

5  ENT1=ENT-ENTZ 
C=ENT1**NN 
IF(M.LT.Z)     GO    TO    1 
IF(M.GT.MAXM)     GO    TO    2 

CP  =  (  ENTHCM-i  )-<3.*ENTH(M)+8.*ENTH(M  +  2)-ENTH(M  + 

3))/<12.*DEL) 
GO    TO    3 

1  C  P  =  I  -  5  0 . *  EN  TH  <  M* 1 ) *9  6  .  *ENTH  (  M  +  2  )  -  7  2 .  *  ENTH  ( iH3  )  + 

2  2.*ENTH(M+4)-6. *ENTH(NI  + 
5))/( 24.*DELJ 
GO    TO    3 

2  CP  =  <6.*ENTH(M-3)-32.*ElMTH(M-2K72.*ENTH(M-l)- 

96 . *E  NTH ( M ) + 50 . *E NTH ( M  + 
1 ) J/(24.*DEL) 

3  FCTN=A*B*C*CP 
RETURN 

END 
DATA 
87  260C.  2000.  401  4 

40.  15C0.  58 

1745.28    -1738.22    -1731.14    -1724.04    -1716.91    -1709.76    - 

17C2.58 
1695.36    -1668.10    -1680*79    -1673.41    -1665.93    -1658.34    - 

1650.61 
1642. 7C   -1634.59    -1626.26    -1617.69    -1603.85    -1599.73    - 

1590.33 
1580.62    -1570.59    -1560.23    -1549.51    -1538.42    -1526.95    - 

1515. G5 
1502.72    -1489.93    -1476.65    -1462.85    -1448.52    -1433.61    - 

1418.11 
1401.96    -1385.18    -1367.71    -1349.51    -1330.56   -1310.83    - 

1 290.29 
1268.90    -1246.64    -1223.43    -1199.38    -1174.32    -1143.27   - 

1121. 20 
1093.09    -1C63.90    -1033.61    -1002.21       -969. b6      -935.94      - 

9C1.04 
-864.93       -827.60 
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